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

      1: #ifndef included_ALE_Numbering_hh
      2: #define included_ALE_Numbering_hh

      4: #ifndef  included_ALE_ParallelMapping_hh
      5: #include <sieve/ParallelMapping.hh>
      6: #endif

      9: namespace ALE {
     10:     // We have a dichotomy between \emph{types}, describing the structure of objects,
     11:     //   and \emph{concepts}, describing the role these objects play in the algorithm.
     12:     //   Below we identify concepts with potential implementing types.
     13:     //
     14:     //   Concept           Type
     15:     //   -------           ----
     16:     //   Overlap           Sifter
     17:     //   Atlas             ConstantSection, UniformSection
     18:     //   Numbering         UniformSection
     19:     //   GlobalOrder       UniformSection
     20:     //
     21:     // We will use factory types to create objects which satisfy a given concept.
     22:   template<typename Point_, typename Value_ = int, typename Alloc_ = malloc_allocator<Point_> >
     23:   class Numbering : public UniformSection<Point_, Value_, 1, Alloc_> {
     24:   public:
     25:     typedef UniformSection<Point_, Value_, 1, Alloc_> base_type;
     26:     typedef typename base_type::point_type point_type;
     27:     typedef typename base_type::value_type value_type;
     28:     typedef typename base_type::atlas_type atlas_type;
     29:   protected:
     30:     int                       _localSize;
     31:     int                      *_offsets;
     32:     std::map<int, point_type> _invOrder;
     33:   public:
     34:     Numbering(MPI_Comm comm, const int debug = 0) : UniformSection<Point_, Value_, 1, Alloc_>(comm, debug), _localSize(0) {
     35:       this->_offsets    = new int[this->commSize()+1];
     36:       this->_offsets[0] = 0;
     37:     };
     38:     virtual ~Numbering() {
     39:       delete [] this->_offsets;
     40:     };
     41:   public: // Sizes
     42:     int        getLocalSize() const {return this->_localSize;};
     43:     void       setLocalSize(const int size) {this->_localSize = size;};
     44:     int        getGlobalSize() const {return this->_offsets[this->commSize()];};
     45:     int        getGlobalOffset(const int p) const {return this->_offsets[p];};
     46:     const int *getGlobalOffsets() const {return this->_offsets;};
     47:     void       setGlobalOffsets(const int offsets[]) {
     48:       for(int p = 0; p <= this->commSize(); ++p) {
     49:         this->_offsets[p] = offsets[p];
     50:       }
     51:     };
     52:   public: // Indices
     53:     virtual int getIndex(const point_type& point) {
     54:       const value_type& idx = this->restrictPoint(point)[0];
     55:       if (idx >= 0) {
     56:         return idx;
     57:       }
     58:       return -(idx+1);
     59:     };
     60:     virtual void setIndex(const point_type& point, const int index) {this->updatePoint(point, &index);};
     61:     virtual bool isLocal(const point_type& point) {return this->restrictPoint(point)[0] >= 0;};
     62:     virtual bool isRemote(const point_type& point) {return this->restrictPoint(point)[0] < 0;};
     63:     point_type getPoint(const int& index) {return this->_invOrder[index];};
     64:     void setPoint(const int& index, const point_type& point) {this->_invOrder[index] = point;};
     65:   };
     66:   template<typename Point_, typename Value_ = ALE::Point>
     67:   class GlobalOrder : public UniformSection<Point_, Value_> {
     68:   public:
     69:     typedef UniformSection<Point_, Value_> base_type;
     70:     typedef typename base_type::point_type point_type;
     71:     typedef typename base_type::value_type value_type;
     72:     typedef typename base_type::atlas_type atlas_type;
     73:   protected:
     74:     int  _localSize;
     75:     int *_offsets;
     76:   public:
     77:     GlobalOrder(MPI_Comm comm, const int debug = 0) : UniformSection<Point_, Value_>(comm, debug), _localSize(0) {
     78:       this->_offsets    = new int[this->commSize()+1];
     79:       this->_offsets[0] = 0;
     80:     };
     81:     ~GlobalOrder() {
     82:       delete [] this->_offsets;
     83:     };
     84:   public: // Sizes
     85:     int        getLocalSize() const {return this->_localSize;};
     86:     void       setLocalSize(const int size) {this->_localSize = size;};
     87:     int        getGlobalSize() const {return this->_offsets[this->commSize()];};
     88:     int        getGlobalOffset(const int p) const {return this->_offsets[p];};
     89:     const int *getGlobalOffsets() const {return this->_offsets;};
     90:     void       setGlobalOffsets(const int offsets[]) {
     91:       for(int p = 0; p <= this->commSize(); ++p) {
     92:         this->_offsets[p] = offsets[p];
     93:       }
     94:     };
     95:   public: // Indices
     96:     virtual int getIndex(const point_type& p) {
     97:       const int idx = this->restrictPoint(p)[0].prefix;
     98:       if (idx >= 0) {
     99:         return idx;
    100:       }
    101:       return -(idx+1);
    102:     };
    103:     virtual void setIndex(const point_type& p, const int index) {
    104:       const value_type idx(index, this->restrictPoint(p)[0].index);
    105:       this->updatePoint(p, &idx);
    106:     };
    107:     virtual bool isLocal(const point_type& p) {return this->restrictPoint(p)[0].prefix >= 0;};
    108:     virtual bool isRemote(const point_type& p) {return this->restrictPoint(p)[0].prefix < 0;};
    109:   };
    110:   template<typename Bundle_, typename Value_ = int, typename Alloc_ = typename Bundle_::alloc_type>
    111:   class NumberingFactory : public ALE::ParallelObject {
    112:   public:
    113:     typedef Bundle_                                         bundle_type;
    114:     typedef Alloc_                                          alloc_type;
    115:     typedef Value_                                          value_type;
    116:     typedef typename bundle_type::sieve_type                sieve_type;
    117:     typedef typename bundle_type::label_type                label_type;
    118:     typedef typename bundle_type::point_type                point_type;
    119:     typedef typename bundle_type::rank_type                 rank_type;
    120:     typedef typename bundle_type::send_overlap_type         send_overlap_type;
    121:     typedef typename bundle_type::recv_overlap_type         recv_overlap_type;
    122:     typedef Numbering<point_type, value_type, alloc_type>   numbering_type;
    123:     typedef typename alloc_type::template rebind<value_type>::other                              value_alloc_type;
    124:     typedef std::map<bundle_type*, std::map<std::string, std::map<int, Obj<numbering_type> > > > numberings_type;
    125:     typedef GlobalOrder<point_type>                         order_type;
    126:     typedef typename order_type::value_type                 oValue_type;
    127:     typedef typename alloc_type::template rebind<oValue_type>::other         oValue_alloc_type;
    128:     typedef std::map<bundle_type*, std::map<std::string, Obj<order_type> > > orders_type;
    129:   protected:
    130:     numberings_type   _localNumberings;
    131:     numberings_type   _numberings;
    132:     orders_type       _localOrders;
    133:     orders_type       _orders;
    134:     orders_type       _ordersBC;
    135:     const value_type  _unknownNumber;
    136:     const oValue_type _unknownOrder;
    137:   protected:
    138:     NumberingFactory(MPI_Comm comm, const int debug = 0) : ALE::ParallelObject(comm, debug), _unknownNumber(-1), _unknownOrder(-1, 0) {};
    139:   public:
    140:     ~NumberingFactory() {};
    141:   public:
    142:     static const Obj<NumberingFactory>& singleton(MPI_Comm comm, const int debug, bool cleanup = false) {
    143:       static Obj<NumberingFactory> *_singleton = NULL;

    145:       if (cleanup) {
    146:         if (debug) {std::cout << "Destroying NumberingFactory" << std::endl;}
    147:         if (_singleton) {delete _singleton;}
    148:         _singleton = NULL;
    149:       } else if (_singleton == NULL) {
    150:         if (debug) {std::cout << "Creating new NumberingFactory" << std::endl;}
    151:         _singleton  = new Obj<NumberingFactory>();
    152:         *_singleton = new NumberingFactory(comm, debug);
    153:       }
    154:       return *_singleton;
    155:     };
    156:     void clear() {
    157:       this->_localNumberings.clear();
    158:       this->_numberings.clear();
    159:       this->_localOrders.clear();
    160:       this->_orders.clear();
    161:       this->_ordersBC.clear();
    162:     };
    163:   public: // Dof ordering
    164:     template<typename Section_>
    165:     void orderPointNew(const Obj<Section_>& section, const Obj<sieve_type>& sieve, const typename Section_::point_type& point, value_type& offset, value_type& bcOffset, const Obj<send_overlap_type>& sendOverlap = NULL) {
    166:       const typename Section_::chart_type& chart = section->getChart();
    167:       int&                                 idx   = section->getIndex(point);

    169:       // If the point does not exist in the chart, throw an error
    170:       if (chart.count(point) == 0) {
    171:         throw ALE::Exception("Unknown point in ordering");
    172:       }
    173:       // If the point has not been ordered
    174:       if (idx == -1) {
    175:         // Recurse to its cover
    176:         const Obj<typename sieve_type::coneSequence>& cone = sieve->cone(point);
    177:         typename sieve_type::coneSequence::iterator   end  = cone->end();

    179:         for(typename sieve_type::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
    180:           if (this->_debug > 1) {std::cout << "    Recursing to " << *c_iter << std::endl;}
    181:           this->orderPoint(section, sieve, *c_iter, offset, bcOffset, sendOverlap);
    182:         }
    183:         const int dim  = section->getFiberDimension(point);
    184:         const int cDim = section->getConstraintDimension(point);
    185:         const int fDim = dim - cDim;

    187:         // If the point has constrained variables
    188:         if (cDim) {
    189:           if (this->_debug > 1) {std::cout << "  Ordering boundary point " << point << " at " << bcOffset << std::endl;}
    190:           section->setIndexBC(point, bcOffset);
    191:           bcOffset += cDim;
    192:         }
    193:         // If the point has free variables
    194:         if (fDim) {
    195:           bool number = true;

    197:           // Maybe use template specialization here
    198:           if (!sendOverlap.isNull() && sendOverlap->capContains(point)) {
    199:             const Obj<typename send_overlap_type::supportSequence>& ranks = sendOverlap->support(point);

    201:             for(typename send_overlap_type::supportSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
    202:               if (this->commRank() > *r_iter) {
    203:                 number = false;
    204:                 break;
    205:               }
    206:             }
    207:           }
    208:           if (number) {
    209:             if (this->_debug > 1) {std::cout << "  Ordering point " << point << " at " << offset << std::endl;}
    210:             section->setIndex(point, offset);
    211:             offset += dim;
    212:           } else {
    213:             if (this->_debug > 1) {std::cout << "  Ignoring ghost point " << point << std::endl;}
    214:           }
    215:         }
    216:       }
    217:     }
    218:     template<typename Section_>
    219:     void orderPoint(const Obj<Section_>& section, const Obj<sieve_type>& sieve, const typename Section_::point_type& point, value_type& offset, value_type& bcOffset, const Obj<send_overlap_type>& sendOverlap = NULL) {
    220:       const Obj<typename Section_::atlas_type>&     atlas = section->getAtlas();
    221:       const Obj<typename sieve_type::coneSequence>& cone = sieve->cone(point);
    222:       typename sieve_type::coneSequence::iterator   end  = cone->end();
    223:       typename Section_::index_type                 idx  = section->getAtlas()->restrictPoint(point)[0];
    224:       const value_type&                             dim  = idx.prefix;
    225:       const typename Section_::index_type           defaultIdx(0, -1);

    227:       if (atlas->getChart().count(point) == 0) {
    228:         idx = defaultIdx;
    229:       }
    230:       if (idx.index == -1) {
    231:         for(typename sieve_type::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
    232:           if (this->_debug > 1) {std::cout << "    Recursing to " << *c_iter << std::endl;}
    233:           this->orderPoint(section, sieve, *c_iter, offset, bcOffset, sendOverlap);
    234:         }
    235:         if (dim > 0) {
    236:           bool number = true;

    238:           // Maybe use template specialization here
    239:           if (!sendOverlap.isNull() && sendOverlap->capContains(point)) {
    240:             const Obj<typename send_overlap_type::supportSequence>& ranks = sendOverlap->support(point);

    242:             for(typename send_overlap_type::supportSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
    243:               if (this->commRank() > *r_iter) {
    244:                 number = false;
    245:                 break;
    246:               }
    247:             }
    248:           }
    249:           if (number) {
    250:             if (this->_debug > 1) {std::cout << "  Ordering point " << point << " at " << offset << std::endl;}
    251:             idx.index = offset;
    252:             atlas->updatePoint(point, &idx);
    253:             offset += dim;
    254:           } else {
    255:             if (this->_debug > 1) {std::cout << "  Ignoring ghost point " << point << std::endl;}
    256:           }
    257:         } else if (dim < 0) {
    258:           if (this->_debug > 1) {std::cout << "  Ordering boundary point " << point << " at " << bcOffset << std::endl;}
    259:           idx.index = bcOffset;
    260:           atlas->updatePoint(point, &idx);
    261:           bcOffset += dim;
    262:         }
    263:       }
    264:     }
    265:     template<typename Section_>
    266:     void orderPatch(const Obj<Section_>& section, const Obj<sieve_type>& sieve, const Obj<send_overlap_type>& sendOverlap = NULL, const value_type offset = 0, const value_type bcOffset = -2) {
    267:       const typename Section_::chart_type& chart = section->getChart();
    268:       int off   = offset;
    269:       int bcOff = bcOffset;

    271:       if (this->_debug > 1) {std::cout << "Ordering patch" << std::endl;}
    272:       for(typename Section_::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
    273:         if (this->_debug > 1) {std::cout << "Ordering closure of point " << *p_iter << std::endl;}
    274:         this->orderPoint(section, sieve, *p_iter, off, bcOff, sendOverlap);
    275:       }
    276:       for(typename Section_::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
    277:         const int& idx  = section->getIndex(*p_iter);

    279:         if (idx < 0) {
    280:           if (this->_debug > 1) {std::cout << "Correcting boundary offset of point " << *p_iter << std::endl;}
    281:           section->setIndex(*p_iter, off - (idx + 2));
    282:         }
    283:       }
    284:     }
    285:   public: // Numbering
    286:     // Number all local points
    287:     //   points in the overlap are only numbered by the owner with the lowest rank
    288:     template<typename Iterator_>
    289:     void constructLocalNumbering(const Obj<numbering_type>& numbering, const Obj<send_overlap_type>& sendOverlap, const Iterator_& pointsBegin, const Iterator_& pointsEnd) {
    290:       const int debug = sendOverlap->debug();
    291:       int localSize = 0;

    293:       if (debug) {std::cout << "["<<numbering->commRank()<<"] Constructing local numbering" << std::endl;}
    294:       for(Iterator_ l_iter = pointsBegin; l_iter != pointsEnd; ++l_iter) {
    295:         numbering->setFiberDimension(*l_iter, 1);
    296:       }
    297:       for(Iterator_ l_iter = pointsBegin; l_iter != pointsEnd; ++l_iter) {
    298:         value_type val;

    300:         if (debug) {std::cout << "["<<numbering->commRank()<<"]   Checking point " << *l_iter << std::endl;}
    301:         if (sendOverlap->capContains(*l_iter)) {
    302:           const Obj<typename send_overlap_type::supportSequence>& sendRanks = sendOverlap->support(*l_iter);
    303:           int minRank = sendOverlap->commSize();

    305:           for(typename send_overlap_type::supportSequence::iterator p_iter = sendRanks->begin(); p_iter != sendRanks->end(); ++p_iter) {
    306:             if (*p_iter < minRank) minRank = *p_iter;
    307:           }
    308:           if (minRank < sendOverlap->commRank()) {
    309:             if (debug) {std::cout << "["<<numbering->commRank()<<"]     remote point, on proc " << minRank << std::endl;}
    310:             val = this->_unknownNumber;
    311:           } else {
    312:             if (debug) {std::cout << "["<<numbering->commRank()<<"]     local point" << std::endl;}
    313:             val = localSize++;
    314:           }
    315:         } else {
    316:           if (debug) {std::cout << "["<<numbering->commRank()<<"]     local point" << std::endl;}
    317:           val = localSize++;
    318:         }
    319:         if (debug) {std::cout << "["<<numbering->commRank()<<"]     has number " << val << std::endl;}
    320:         numbering->updatePoint(*l_iter, &val);
    321:       }
    322:       if (debug) {std::cout << "["<<numbering->commRank()<<"]   local points" << std::endl;}
    323:       numbering->setLocalSize(localSize);
    324:     }
    325:     // Order all local points
    326:     //   points in the overlap are only ordered by the owner with the lowest rank
    327:     template<typename Iterator_, typename Section_>
    328:     void constructLocalOrder(const Obj<order_type>& order, const Obj<send_overlap_type>& sendOverlap, const Iterator_& pointsBegin, const Iterator_& pointsEnd, const Obj<Section_>& section, const int space = -1, const bool withBC = false, const Obj<label_type>& label = NULL) {
    329:       const int debug = sendOverlap->debug();
    330:       int localSize = 0;

    332:       if (debug) {std::cout << "["<<order->commRank()<<"] Constructing local ordering" << std::endl;}
    333:       for(Iterator_ l_iter = pointsBegin; l_iter != pointsEnd; ++l_iter) {
    334:         order->setFiberDimension(*l_iter, 1);
    335:       }
    336:       for(Iterator_ l_iter = pointsBegin; l_iter != pointsEnd; ++l_iter) {
    337:         oValue_type val;

    339:         if (debug) {std::cout << "["<<order->commRank()<<"]   Checking point " << *l_iter << std::endl;}
    340:         if (sendOverlap->capContains(*l_iter)) {
    341:           const Obj<typename send_overlap_type::supportSequence>& sendPatches = sendOverlap->support(*l_iter);
    342:           int minRank = sendOverlap->commSize();

    344:           for(typename send_overlap_type::supportSequence::iterator p_iter = sendPatches->begin(); p_iter != sendPatches->end(); ++p_iter) {
    345:             if (*p_iter < minRank) minRank = *p_iter;
    346:           }
    347:           bool remotePoint = (minRank < sendOverlap->commRank()) || (!label.isNull() && (label->cone(*l_iter)->size() > 0));

    349:           if (remotePoint) {
    350:             if (debug) {std::cout << "["<<order->commRank()<<"]     remote point, on proc " << minRank << std::endl;}
    351:             val = this->_unknownOrder;
    352:           } else {
    353:             if (debug) {std::cout << "["<<order->commRank()<<"]     local point" << std::endl;}
    354:             val.prefix = localSize;
    355:             if (withBC) {
    356:               val.index  = space < 0 ? section->getFiberDimension(*l_iter) : section->getFiberDimension(*l_iter, space);
    357:             } else {
    358:               val.index  = space < 0 ? section->getConstrainedFiberDimension(*l_iter) : section->getConstrainedFiberDimension(*l_iter, space);
    359:             }
    360:           }
    361:         } else {
    362:           if (debug) {std::cout << "["<<order->commRank()<<"]     local point" << std::endl;}
    363:           val.prefix = localSize;
    364:           if (withBC) {
    365:             val.index  = space < 0 ? section->getFiberDimension(*l_iter) : section->getFiberDimension(*l_iter, space);
    366:           } else {
    367:             val.index  = space < 0 ? section->getConstrainedFiberDimension(*l_iter) : section->getConstrainedFiberDimension(*l_iter, space);
    368:           }
    369:         }
    370:         if (debug) {std::cout << "["<<order->commRank()<<"]     has offset " << val.prefix << " and size " << val.index << std::endl;}
    371:         localSize += val.index;
    372:         order->updatePoint(*l_iter, &val);
    373:       }
    374:       if (debug) {std::cout << "["<<order->commRank()<<"]   local size" << localSize << std::endl;}
    375:       order->setLocalSize(localSize);
    376:     }
    377:     // Calculate process offsets
    378:     template<typename Numbering>
    379:     void calculateOffsets(const Obj<Numbering>& numbering) {
    380:       int  localSize = numbering->getLocalSize();
    381:       int *offsets   = new int[numbering->commSize()+1];

    383:       offsets[0] = 0;
    384:       MPI_Allgather(&localSize, 1, MPI_INT, &(offsets[1]), 1, MPI_INT, numbering->comm());
    385:       for(int p = 2; p <= numbering->commSize(); p++) {
    386:         offsets[p] += offsets[p-1];
    387:       }
    388:       numbering->setGlobalOffsets(offsets);
    389:       delete [] offsets;
    390:     }
    391:     // Update local offsets based upon process offsets
    392:     template<typename Numbering, typename Iterator>
    393:     void updateOrder(const Obj<Numbering>& numbering, const Iterator& pointsBegin, const Iterator& pointsEnd) {
    394:       const typename Numbering::value_type val = numbering->getGlobalOffset(numbering->commRank());

    396:       for(Iterator l_iter = pointsBegin; l_iter != pointsEnd; ++l_iter) {
    397:         if (numbering->isLocal(*l_iter)) {
    398:           numbering->updateAddPoint(*l_iter, &val);
    399:         }
    400:       }
    401:     }
    402:     // Communicate numbers in the overlap
    403:     void completeNumbering(const Obj<numbering_type>& numbering, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, bool allowDuplicates = false) {
    404: #if 1
    405:       typedef ALE::UniformSection<ALE::Pair<int, typename send_overlap_type::source_type>, typename numbering_type::value_type> OverlapSection;
    406:       typedef typename OverlapSection::point_type overlap_point_type;
    407:       Obj<OverlapSection> overlapSection = new OverlapSection(numbering->comm(), numbering->debug());
    408:       const int debug = sendOverlap->debug();

    410:       if (debug) {std::cout << "["<<numbering->commRank()<<"] Completing numbering" << std::endl;}
    411:       ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, numbering, overlapSection);
    412:       if (debug) {overlapSection->view("Overlap Section");}
    413:       const typename recv_overlap_type::capSequence::iterator rBegin = recvOverlap->capBegin();
    414:       const typename recv_overlap_type::capSequence::iterator rEnd   = recvOverlap->capEnd();

    416:       for(typename recv_overlap_type::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
    417:         const int                                             rank    = *r_iter;
    418:         const typename recv_overlap_type::supportSequence::iterator pBegin  = recvOverlap->supportBegin(*r_iter);
    419:         const typename recv_overlap_type::supportSequence::iterator pEnd    = recvOverlap->supportEnd(*r_iter);

    421:         for(typename recv_overlap_type::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
    422:           const point_type& localPoint  = *p_iter;
    423:           const point_type& remotePoint = p_iter.color();
    424:           const overlap_point_type                   oPoint = overlap_point_type(rank, remotePoint);
    425:           const int                                  fDim   = overlapSection->getFiberDimension(oPoint);
    426:           const typename numbering_type::value_type *values = overlapSection->restrictPoint(oPoint);

    428:           for(int i = 0; i < fDim; ++i) {
    429:             if (debug) {std::cout << "["<<numbering->commRank()<<"]     local point " << localPoint << " remote point " << remotePoint << " number " << values[i] << std::endl;}
    430:             if (values[i] >= 0) {
    431:               if (numbering->isLocal(localPoint) && !allowDuplicates) {
    432:                 ostringstream msg;
    433:                 msg << "["<<numbering->commRank()<<"]Multiple indices for local point " << localPoint << " remote point " << remotePoint << " from " << rank << " with index " << values[i];
    434:                 throw ALE::Exception(msg.str().c_str());
    435:               }
    436:               if (!numbering->hasPoint(localPoint)) {
    437:                 ostringstream msg;
    438:                 msg << "["<<numbering->commRank()<<"]Unexpected local point " << localPoint << " remote point " << remotePoint << " from " << rank << " with index " << values[i];
    439:                 throw ALE::Exception(msg.str().c_str());
    440:               }
    441:               int val = -(values[i]+1);
    442:               numbering->updatePoint(localPoint, &val);
    443:             }
    444:           }
    445:         }
    446:       }
    447: #else
    448:       typedef Field<send_overlap_type, int, Section<point_type, value_type, value_alloc_type> > send_section_type;
    449:       typedef Field<recv_overlap_type, int, Section<point_type, value_type, value_alloc_type> > recv_section_type;
    450:       typedef typename ALE::DiscreteSieve<point_type, alloc_type>                   dsieve_type;
    451:       typedef typename ALE::Topology<int, dsieve_type, alloc_type>                  dtopology_type;
    452:       typedef typename ALE::New::SectionCompletion<dtopology_type, int, alloc_type> completion;
    453:       const Obj<send_section_type> sendSection = new send_section_type(numbering->comm(), this->debug());
    454:       const Obj<recv_section_type> recvSection = new recv_section_type(numbering->comm(), sendSection->getTag(), this->debug());
    455:       const int debug = sendOverlap->debug();

    457:       if (debug) {std::cout << "["<<numbering->commRank()<<"] Completing numbering" << std::endl;}
    458:       completion::completeSection(sendOverlap, recvOverlap, numbering->getAtlas(), numbering, sendSection, recvSection);
    459:       const Obj<typename recv_overlap_type::baseSequence> rPoints = recvOverlap->base();

    461:       for(typename recv_overlap_type::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rPoints->end(); ++p_iter) {
    462:         const Obj<typename recv_overlap_type::coneSequence>& ranks      = recvOverlap->cone(*p_iter);
    463:         const typename recv_overlap_type::target_type&       localPoint = *p_iter;

    465:         for(typename recv_overlap_type::coneSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
    466:           const typename recv_overlap_type::target_type&       remotePoint = r_iter.color();
    467:           const int                                            rank        = *r_iter;
    468:           const Obj<typename recv_section_type::section_type>& section     = recvSection->getSection(rank);
    469:           const typename recv_section_type::value_type        *values      = section->restrictPoint(remotePoint);

    471:           if (section->getFiberDimension(remotePoint) == 0) continue;
    472:           if (debug) {std::cout << "["<<numbering->commRank()<<"]     local point " << localPoint << " remote point " << remotePoint << " number " << values[0] << std::endl;}
    473:           if (values[0] >= 0) {
    474:             if (debug) {std::cout << "["<<numbering->commRank()<<"] local point " << localPoint << " dim " << numbering->getAtlas()->getFiberDimension(localPoint) << std::endl;}
    475:             if (numbering->isLocal(localPoint) && !allowDuplicates) {
    476:               ostringstream msg;
    477:               msg << "["<<numbering->commRank()<<"]Multiple indices for local point " << localPoint << " remote point " << remotePoint << " from " << rank << " with index " << values[0];
    478:               throw ALE::Exception(msg.str().c_str());
    479:             }
    480:             if (numbering->getAtlas()->getFiberDimension(localPoint) == 0) {
    481:               ostringstream msg;
    482:               msg << "["<<numbering->commRank()<<"]Unexpected local point " << localPoint << " remote point " << remotePoint << " from " << rank << " with index " << values[0];
    483:               throw ALE::Exception(msg.str().c_str());
    484:             }
    485:             int val = -(values[0]+1);
    486:             numbering->updatePoint(localPoint, &val);
    487:           }
    488:         }
    489:       }
    490: #endif
    491:     }
    492:     // Communicate (size,offset)s in the overlap
    493:     void completeOrder(const Obj<order_type>& order, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, bool allowDuplicates = false) {
    494: #if 1
    495:       typedef ALE::UniformSection<ALE::Pair<int, typename send_overlap_type::source_type>, typename order_type::value_type> OverlapSection;
    496:       typedef typename OverlapSection::point_type overlap_point_type;
    497:       Obj<OverlapSection> overlapSection = new OverlapSection(order->comm(), order->debug());
    498:       const int debug = sendOverlap->debug();

    500:       if (debug) {std::cout << "["<<order->commRank()<<"] Completing ordering" << std::endl;}
    501:       ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, order, overlapSection);
    502:       if (debug) {overlapSection->view("Overlap Section");}
    503:       const typename recv_overlap_type::capSequence::iterator rBegin = recvOverlap->capBegin();
    504:       const typename recv_overlap_type::capSequence::iterator rEnd   = recvOverlap->capEnd();

    506:       for(typename recv_overlap_type::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
    507:         const int                                             rank    = *r_iter;
    508:         const typename recv_overlap_type::supportSequence::iterator pBegin  = recvOverlap->supportBegin(*r_iter);
    509:         const typename recv_overlap_type::supportSequence::iterator pEnd    = recvOverlap->supportEnd(*r_iter);

    511:         for(typename recv_overlap_type::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
    512:           const point_type& localPoint  = *p_iter;
    513:           const point_type& remotePoint = p_iter.color();
    514:           const overlap_point_type               oPoint = overlap_point_type(rank, remotePoint);
    515:           const int                              fDim   = overlapSection->getFiberDimension(oPoint);
    516:           const typename order_type::value_type *values = overlapSection->restrictPoint(oPoint);

    518:           for(int i = 0; i < fDim; ++i) {
    519:             if (debug) {std::cout << "["<<order->commRank()<<"]     local point " << localPoint << " remote point " << remotePoint<<"("<<rank<<")" << " offset " << values[i].prefix << " and size " << values[i].index << std::endl;}
    520:             if (values[i].index == 0) continue;
    521:             if (values[i].prefix >= 0) {
    522:               if (order->isLocal(localPoint)) {
    523:                 if (!allowDuplicates) {
    524:                   ostringstream msg;
    525:                   msg << "["<<order->commRank()<<"]Multiple indices for local point " << localPoint << " remote point " << remotePoint << " from " << rank << " with index " << values[i];
    526:                   throw ALE::Exception(msg.str().c_str());
    527:                 }
    528:                 continue;
    529:               }
    530:               const oValue_type val(-(values[i].prefix+1), values[i].index);
    531:               order->updatePoint(localPoint, &val);
    532:             } else {
    533:               if (order->isLocal(localPoint)) continue;
    534:               order->updatePoint(localPoint, &values[i]);
    535:             }
    536:           }
    537:         }
    538:       }
    539: #else
    540:       typedef Field<send_overlap_type, int, Section<point_type, oValue_type, oValue_alloc_type> > send_section_type;
    541:       typedef Field<recv_overlap_type, int, Section<point_type, oValue_type, oValue_alloc_type> > recv_section_type;
    542:       typedef ConstantSection<point_type, int, alloc_type>                          constant_sizer;
    543:       typedef typename ALE::DiscreteSieve<point_type, alloc_type>                   dsieve_type;
    544:       typedef typename ALE::Topology<int, dsieve_type, alloc_type>                  dtopology_type;
    545:       typedef typename ALE::New::SectionCompletion<dtopology_type, int, alloc_type> completion;
    546:       const Obj<send_section_type> sendSection = new send_section_type(order->comm(), this->debug());
    547:       const Obj<recv_section_type> recvSection = new recv_section_type(order->comm(), sendSection->getTag(), this->debug());
    548:       const int debug = sendOverlap->debug();

    550:       if (debug) {std::cout << "["<<order->commRank()<<"] Completing ordering" << std::endl;}
    551:       completion::completeSection(sendOverlap, recvOverlap, order->getAtlas(), order, sendSection, recvSection);
    552:       Obj<typename recv_overlap_type::baseSequence> recvPoints = recvOverlap->base();

    554:       for(typename recv_overlap_type::baseSequence::iterator p_iter = recvPoints->begin(); p_iter != recvPoints->end(); ++p_iter) {
    555:         if (!order->hasPoint(*p_iter)) {
    556:           order->setFiberDimension(*p_iter, 1);
    557:           order->updatePoint(*p_iter, &this->_unknownOrder);
    558:         }
    559:       }
    560:       for(typename recv_overlap_type::baseSequence::iterator p_iter = recvPoints->begin(); p_iter != recvPoints->end(); ++p_iter) {
    561:         const Obj<typename recv_overlap_type::coneSequence>& ranks      = recvOverlap->cone(*p_iter);
    562:         const typename recv_overlap_type::target_type&       localPoint = *p_iter;

    564:         for(typename recv_overlap_type::coneSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
    565:           const typename recv_overlap_type::target_type&       remotePoint = r_iter.color();
    566:           const int                                            rank        = *r_iter;
    567:           const Obj<typename recv_section_type::section_type>& section     = recvSection->getSection(rank);
    568:           const typename recv_section_type::value_type        *values      = section->restrictPoint(remotePoint);

    570:           if (section->getFiberDimension(remotePoint) == 0) continue;
    571:           if (debug) {std::cout << "["<<order->commRank()<<"]     local point " << localPoint << " remote point " << remotePoint<<"("<<rank<<")" << " offset " << values[0].prefix << " and size " << values[0].index << std::endl;}
    572:           if (values[0].index == 0) continue;
    573:           if (values[0].prefix >= 0) {
    574:             if (order->isLocal(localPoint)) {
    575:               if (!allowDuplicates) {
    576:                 ostringstream msg;
    577:                 msg << "["<<order->commRank()<<"]Multiple indices for local point " << localPoint << " remote point " << remotePoint << " from " << rank << " with index " << values[0];
    578:                 throw ALE::Exception(msg.str().c_str());
    579:               }
    580:               continue;
    581:             }
    582:             const oValue_type val(-(values[0].prefix+1), values[0].index);
    583:             order->updatePoint(localPoint, &val);
    584:           } else {
    585:             if (order->isLocal(localPoint)) continue;
    586:             order->updatePoint(localPoint, values);
    587:           }
    588:         }
    589:       }
    590: #endif
    591:     }
    592:     // Construct a full global numbering
    593:     template<typename Iterator>
    594:     void constructNumbering(const Obj<numbering_type>& numbering, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Iterator& pointsBegin, const Iterator& pointsEnd) {
    595:       this->constructLocalNumbering(numbering, sendOverlap, pointsBegin, pointsEnd);
    596:       this->calculateOffsets(numbering);
    597:       this->updateOrder(numbering, pointsBegin, pointsEnd);
    598:       this->completeNumbering(numbering, sendOverlap, recvOverlap);
    599:     }
    600:     // Construct a full global order
    601:     template<typename Iterator, typename Section>
    602:     void constructOrder(const Obj<order_type>& order, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Iterator& pointsBegin, const Iterator& pointsEnd, const Obj<Section>& section, const int space = -1, const Obj<label_type>& label = NULL) {
    603:       this->constructLocalOrder(order, sendOverlap, pointsBegin, pointsEnd, section, space, false, label);
    604:       this->calculateOffsets(order);
    605:       this->updateOrder(order, pointsBegin, pointsEnd);
    606:       this->completeOrder(order, sendOverlap, recvOverlap);
    607:     }
    608:     template<typename Iterator, typename Section>
    609:     void constructOrderBC(const Obj<order_type>& order, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Iterator& pointsBegin, const Iterator& pointsEnd, const Obj<Section>& section, const int space = -1, const Obj<label_type>& label = NULL) {
    610:       this->constructLocalOrder(order, sendOverlap, pointsBegin, pointsEnd, section, space, true, label);
    611:       this->calculateOffsets(order);
    612:       this->updateOrder(order, pointsBegin, pointsEnd);
    613:       this->completeOrder(order, sendOverlap, recvOverlap);
    614:     }
    615:   public:
    616:     // Construct the inverse map from numbers to points
    617:     //   If we really need this, then we should consider using a label
    618:     void constructInverseOrder(const Obj<numbering_type>& numbering) {
    619:       const typename numbering_type::chart_type& chart = numbering->getChart();

    621:       for(typename numbering_type::chart_type::iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
    622:         numbering->setPoint(numbering->getIndex(*p_iter), *p_iter);
    623:       }
    624:     }
    625:   public: // Real interface
    626:     template<typename ABundle_>
    627:     const Obj<numbering_type>& getLocalNumbering(const Obj<ABundle_>& bundle, const int depth) {
    628:       if ((this->_localNumberings.find(bundle.ptr()) == this->_localNumberings.end()) ||
    629:           (this->_localNumberings[bundle.ptr()].find("depth") == this->_localNumberings[bundle.ptr()].end()) ||
    630:           (this->_localNumberings[bundle.ptr()]["depth"].find(depth) == this->_localNumberings[bundle.ptr()]["depth"].end())) {
    631:         Obj<numbering_type>    numbering   = new numbering_type(bundle->comm(), bundle->debug());
    632:         Obj<send_overlap_type> sendOverlap = new send_overlap_type(bundle->comm(), bundle->debug());

    634:         this->constructLocalNumbering(numbering, sendOverlap, bundle->depthStratum(depth)->begin(), bundle->depthStratum(depth)->end());
    635:         if (this->_debug) {std::cout << "Creating new local numbering: ptr " << bundle.ptr() << " depth " << depth << std::endl;}
    636:         this->_localNumberings[bundle.ptr()]["depth"][depth] = numbering;
    637:       } else {
    638:         if (this->_debug) {std::cout << "Using old local numbering: ptr " << bundle.ptr() << " depth " << depth << std::endl;}
    639:       }
    640:       return this->_localNumberings[bundle.ptr()]["depth"][depth];
    641:     }
    642:     template<typename ABundle_>
    643:     const Obj<numbering_type>& getNumbering(const Obj<ABundle_>& bundle, const int depth) {
    644:       if ((this->_numberings.find(bundle.ptr()) == this->_numberings.end()) ||
    645:           (this->_numberings[bundle.ptr()].find("depth") == this->_numberings[bundle.ptr()].end()) ||
    646:           (this->_numberings[bundle.ptr()]["depth"].find(depth) == this->_numberings[bundle.ptr()]["depth"].end())) {
    647:         bundle->constructOverlap();
    648:         Obj<numbering_type>    numbering   = new numbering_type(bundle->comm(), bundle->debug());
    649:         Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
    650:         Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();

    652:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Creating new numbering: fixed depth value " << depth << std::endl;}
    653:         if (depth == -1) {
    654:           this->constructNumbering(numbering, sendOverlap, recvOverlap, bundle->getSieve()->getChart().begin(), bundle->getSieve()->getChart().end());
    655:         } else {
    656:           this->constructNumbering(numbering, sendOverlap, recvOverlap, bundle->depthStratum(depth)->begin(), bundle->depthStratum(depth)->end());
    657:         }
    658:         this->_numberings[bundle.ptr()]["depth"][depth] = numbering;
    659:       } else {
    660:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Using old numbering: fixed depth value " << depth << std::endl;}
    661:       }
    662:       return this->_numberings[bundle.ptr()]["depth"][depth];
    663:     }
    664:     template<typename ABundle_>
    665:     const Obj<numbering_type>& getNumbering(const Obj<ABundle_>& bundle, const std::string& labelname, const int value) {
    666:       if ((this->_numberings.find(bundle.ptr()) == this->_numberings.end()) ||
    667:           (this->_numberings[bundle.ptr()].find(labelname) == this->_numberings[bundle.ptr()].end()) ||
    668:           (this->_numberings[bundle.ptr()][labelname].find(value) == this->_numberings[bundle.ptr()][labelname].end())) {
    669:         bundle->constructOverlap();
    670:         Obj<numbering_type>    numbering   = new numbering_type(bundle->comm(), bundle->debug());
    671:         Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
    672:         Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();

    674:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Creating new numbering: " << labelname << " value " << value << std::endl;}
    675:         this->constructNumbering(numbering, sendOverlap, recvOverlap, bundle->getLabelStratum(labelname, value)->begin(), bundle->getLabelStratum(labelname, value)->end());
    676:         this->_numberings[bundle.ptr()][labelname][value] = numbering;
    677:       } else {
    678:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Using old numbering: " << labelname << " value " << value << std::endl;}
    679:       }
    680:       return this->_numberings[bundle.ptr()][labelname][value];
    681:     }
    682:     template<typename ABundle_, typename Section_>
    683:     const Obj<order_type>& getLocalOrder(const Obj<ABundle_>& bundle, const std::string& name, const Obj<Section_>& section, const int space = -1, const Obj<label_type>& label = NULL) {
    684:       if ((this->_localOrders.find(bundle.ptr()) == this->_localOrders.end()) ||
    685:           (this->_localOrders[bundle.ptr()].find(name) == this->_localOrders[bundle.ptr()].end())) {
    686:         Obj<order_type>        order       = new order_type(bundle->comm(), bundle->debug());
    687:         Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
    688:         Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();

    690:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Creating new local order: " << name << std::endl;}
    691:         this->constructLocalOrder(order, sendOverlap, section->getChart().begin(), section->getChart().end(), section, space, false, label);
    692:         this->_localOrders[bundle.ptr()][name] = order;
    693:       } else {
    694:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Using old local order: " << name << std::endl;}
    695:       }
    696:       return this->_localOrders[bundle.ptr()][name];
    697:     }
    698:     template<typename ABundle_, typename Section_>
    699:     const Obj<order_type>& getGlobalOrder(const Obj<ABundle_>& bundle, const std::string& name, const Obj<Section_>& section, const int space = -1, const Obj<label_type>& label = NULL) {
    700:       if ((this->_orders.find(bundle.ptr()) == this->_orders.end()) ||
    701:           (this->_orders[bundle.ptr()].find(name) == this->_orders[bundle.ptr()].end())) {
    702:         bundle->constructOverlap();
    703:         Obj<order_type>        order       = new order_type(bundle->comm(), bundle->debug());
    704:         Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
    705:         Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();

    707:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Creating new global order: " << name << std::endl;}
    708:         this->constructOrder(order, sendOverlap, recvOverlap, section->getChart().begin(), section->getChart().end(), section, space, label);
    709:         this->_orders[bundle.ptr()][name] = order;
    710:       } else {
    711:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Using old global order: " << name << std::endl;}
    712:       }
    713:       return this->_orders[bundle.ptr()][name];
    714:     }
    715:     template<typename ABundle_, typename Iterator_, typename Section_>
    716:     const Obj<order_type>& getGlobalOrder(const Obj<ABundle_>& bundle, const std::string& name, const Iterator_& pointsBegin, const Iterator_& pointsEnd, const Obj<Section_>& section, const int space = -1, const Obj<label_type>& label = NULL) {
    717:       if ((this->_orders.find(bundle.ptr()) == this->_orders.end()) ||
    718:           (this->_orders[bundle.ptr()].find(name) == this->_orders[bundle.ptr()].end())) {
    719:         bundle->constructOverlap();
    720:         Obj<order_type>        order       = new order_type(bundle->comm(), bundle->debug());
    721:         Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
    722:         Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();

    724:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Creating new global order: " << name << std::endl;}
    725:         this->constructOrder(order, sendOverlap, recvOverlap, pointsBegin, pointsEnd, section, space, label);
    726:         this->_orders[bundle.ptr()][name] = order;
    727:       } else {
    728:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Using old global order: " << name << std::endl;}
    729:       }
    730:       return this->_orders[bundle.ptr()][name];
    731:     }
    732:     template<typename ABundle_, typename Section_>
    733:     const Obj<order_type>& getGlobalOrderWithBC(const Obj<ABundle_>& bundle, const std::string& name, const Obj<Section_>& section, const int space = -1, const Obj<label_type>& label = NULL) {
    734:       if ((this->_orders.find(bundle.ptr()) == this->_orders.end()) ||
    735:           (this->_orders[bundle.ptr()].find(name) == this->_orders[bundle.ptr()].end())) {
    736:         bundle->constructOverlap();
    737:         Obj<order_type>        order       = new order_type(bundle->comm(), bundle->debug());
    738:         Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
    739:         Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();

    741:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Creating new global order: " << name << std::endl;}
    742:         this->constructOrderBC(order, sendOverlap, recvOverlap, section->getChart().begin(), section->getChart().end(), section, space, label);
    743:         this->_orders[bundle.ptr()][name] = order;
    744:       } else {
    745:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Using old global order: " << name << std::endl;}
    746:       }
    747:       return this->_orders[bundle.ptr()][name];
    748:     }
    749:     template<typename ABundle_, typename Iterator_, typename Section_>
    750:     const Obj<order_type>& getGlobalOrderWithBC(const Obj<ABundle_>& bundle, const std::string& name, const Iterator_& pointsBegin, const Iterator_& pointsEnd, const Obj<Section_>& section, const int space = -1, const Obj<label_type>& label = NULL) {
    751:       if ((this->_ordersBC.find(bundle.ptr()) == this->_ordersBC.end()) ||
    752:           (this->_ordersBC[bundle.ptr()].find(name) == this->_ordersBC[bundle.ptr()].end())) {
    753:         bundle->constructOverlap();
    754:         Obj<order_type>        order       = new order_type(bundle->comm(), bundle->debug());
    755:         Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
    756:         Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();

    758:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Creating new global order: " << name << std::endl;}
    759:         this->constructOrderBC(order, sendOverlap, recvOverlap, pointsBegin, pointsEnd, section, space, label);
    760:         this->_orders[bundle.ptr()][name] = order;
    761:       } else {
    762:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Using old global order: " << name << std::endl;}
    763:       }
    764:       return this->_orders[bundle.ptr()][name];
    765:     }
    766:     template<typename ABundle_>
    767:     void setGlobalOrder(const Obj<ABundle_>& bundle, const std::string& name, const Obj<order_type>& order) {
    768:       this->_orders[bundle.ptr()][name] = order;
    769:     }
    770:   };
    771: }
    772: #endif