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

      1: #ifndef included_ALE_ParallelMapping_hh
      2: #define included_ALE_ParallelMapping_hh

      4: #ifndef  included_ALE_BasicCommunication_hh
      5: #include <sieve/BasicCommunication.hh>
      6: #endif

      8: #ifndef  included_ALE_IField_hh
      9: #include <sieve/IField.hh>
     10: #endif

     12: #ifndef  included_ALE_Sections_hh
     13: #include <sieve/Sections.hh>
     14: #endif

     16: #include <functional>
     17: #include <valarray>

     19: namespace ALE {
     20:   template<class _Tp>
     21:   struct Identity : public std::unary_function<_Tp,_Tp>
     22:   {
     23:     _Tp& operator()(_Tp& x) const {return x;}
     24:     const _Tp& operator()(const _Tp& x) const {return x;}
     25:   };

     27:   template<class _Tp>
     28:   struct IsEqual : public std::unary_function<_Tp, bool>, public std::binary_function<_Tp, _Tp, bool>
     29:   {
     30:     const _Tp& x;
     31:     IsEqual(const _Tp& x) : x(x) {};
     32:     bool operator()(_Tp& y) const {return x == y;}
     33:     bool operator()(const _Tp& y) const {return x == y;}
     34:     bool operator()(_Tp& y, _Tp& dummy) const {return x == y;}
     35:     bool operator()(const _Tp& y, const _Tp& dummy) const {return x == y;}
     36:   };

     38:   // Creates new global point names and renames local points globally
     39:   template<typename Point>
     40:   class PointFactory : ALE::ParallelObject {
     41:   public:
     42:     typedef Point                           point_type;
     43:     typedef std::map<point_type,point_type> renumbering_type;
     44:     typedef std::map<int,std::map<point_type,point_type> > remote_renumbering_type;
     45:   protected:
     46:     point_type       originalMax;
     47:     point_type       currentMax;
     48:     renumbering_type renumbering;
     49:     renumbering_type invRenumbering;
     50:     remote_renumbering_type remoteRenumbering;
     51:   protected:
     52:     PointFactory(MPI_Comm comm, const int debug = 0) : ALE::ParallelObject(comm, debug), originalMax(-1) {};
     53:   public:
     54:     ~PointFactory() {};
     55:   public:
     56:     static PointFactory& singleton(MPI_Comm comm, const point_type& maxPoint, const int debug = 0, bool cleanup = false) {
     57:       static PointFactory *_singleton = NULL;

     59:       if (cleanup) {
     60:         if (debug) {std::cout << "Destroying PointFactory" << std::endl;}
     61:         if (_singleton) {delete _singleton;}
     62:         _singleton = NULL;
     63:       } else if (_singleton == NULL) {
     64:         if (debug) {std::cout << "Creating new PointFactory" << std::endl;}
     65:         _singleton  = new PointFactory(comm, debug);
     66:         _singleton->setMax(maxPoint);
     67:       }
     68:       return *_singleton;
     69:     };
     70:     void setMax(const point_type& maxPoint) {
     71:       PetscErrorCode MPI_Allreduce((void *) &maxPoint, &this->originalMax, 1, MPI_INT, MPI_MAX, this->comm());CHKERRXX(ierr);
     72:       ++this->originalMax;
     73:       this->currentMax = this->originalMax;
     74:     };
     75:     void clear() {
     76:       this->currentMax = this->originalMax;
     77:       this->renumbering.clear();
     78:       this->invRenumbering.clear();
     79:     };
     80:   public:
     81:     template<typename Iterator>
     82:     void renumberPoints(const Iterator& begin, const Iterator& end) {
     83:       renumberPoints(begin, end, Identity<typename Iterator::value_type>());
     84:     }
     85:     template<typename Iterator, typename KeyExtractor>
     86:     void renumberPoints(const Iterator& begin, const Iterator& end, const KeyExtractor& ex) {
     87:       int numPoints = 0, numGlobalPoints, firstPoint;

     89:       for(Iterator p_iter = begin; p_iter != end; ++p_iter) ++numPoints;
     90:       MPI_Allreduce(&numPoints, &numGlobalPoints, 1, MPI_INT, MPI_SUM, this->comm());
     91:       MPI_Scan(&numPoints, &firstPoint, 1, MPI_INT, MPI_SUM, this->comm());
     92:       firstPoint += this->currentMax - numPoints;
     93:       this->currentMax += numGlobalPoints;
     94:       for(Iterator p_iter = begin; p_iter != end; ++p_iter, ++firstPoint) {
     95:         if (this->debug()) {std::cout << "["<<this->commRank()<<"]: New point " << ex(*p_iter) << " --> " << firstPoint << std::endl;}
     96:         this->renumbering[firstPoint]     = ex(*p_iter);
     97:         this->invRenumbering[ex(*p_iter)] = firstPoint;
     98:       }
     99:     }
    100:   public:
    101:     void constructRemoteRenumbering() {
    102:       const int localSize   = this->renumbering.size();
    103:       int      *remoteSizes = new int[this->commSize()];
    104:       int      *localMap    = new int[localSize*2];
    105:       int      *recvCounts  = new int[this->commSize()];
    106:       int      *displs      = new int[this->commSize()];

    108:       // Populate arrays
    109:       int r = 0;
    110:       for(typename renumbering_type::const_iterator r_iter = renumbering.begin(); r_iter != renumbering.end(); ++r_iter, ++r) {
    111:         localMap[r*2+0] = r_iter->first;
    112:         localMap[r*2+1] = r_iter->second;
    113:       }
    114:       // Communicate renumbering sizes
    115:       MPI_Allgather((void*) &localSize, 1, MPI_INT, remoteSizes, 1, MPI_INT, this->comm());
    116:       for(int p = 0; p < this->commSize(); ++p) {
    117:         recvCounts[p] = remoteSizes[p]*2;
    118:         if (p == 0) {
    119:           displs[p]   = 0;
    120:         } else {
    121:           displs[p]   = displs[p-1] + recvCounts[p-1];
    122:         }
    123:       }
    124:       int *remoteMaps = new int[displs[this->commSize()-1]+recvCounts[this->commSize()-1]];
    125:       // Communicate renumberings
    126:       MPI_Allgatherv(localMap, localSize*2, MPI_INT, remoteMaps, recvCounts, displs, MPI_INT, this->comm());
    127:       // Populate maps
    128:       for(int p = 0; p < this->commSize(); ++p) {
    129:         if (p == this->commRank()) continue;
    130:         int offset = displs[p];

    132:         for(int r = 0; r < remoteSizes[p]; ++r) {
    133:           this->remoteRenumbering[p][remoteMaps[r*2+0+offset]] = remoteMaps[r*2+1+offset];
    134:           if (this->debug()) {std::cout << "["<<this->commRank()<<"]: Remote renumbering["<<p<<"] " << remoteMaps[r*2+0+offset] << " --> " << remoteMaps[r*2+1+offset] << std::endl;}
    135:         }
    136:       }
    137:       // Cleanup
    138:       delete [] recvCounts;
    139:       delete [] displs;
    140:       delete [] localMap;
    141:       delete [] remoteMaps;
    142:       delete [] remoteSizes;
    143:     };
    144:   public:
    145:     // global point --> local point
    146:     renumbering_type& getRenumbering() {
    147:       return this->renumbering;
    148:     };
    149:     // local point --> global point
    150:     renumbering_type& getInvRenumbering() {
    151:       return this->invRenumbering;
    152:     };
    153:     // rank --> global point --> local point
    154:     remote_renumbering_type& getRemoteRenumbering() {
    155:       return this->remoteRenumbering;
    156:     };
    157:   };

    159:   template<typename Alloc_ = malloc_allocator<int> >
    160:   class OverlapBuilder {
    161:   public:
    162:     typedef Alloc_ alloc_type;
    163:   protected:
    164:     template<typename T>
    165:     struct lessPair: public std::binary_function<std::pair<T,T>, std::pair<T,T>, bool> {
    166:       bool operator()(const std::pair<T,T>& x, const std::pair<T,T>& y) const {
    167:         return x.first < y.first;
    168:       }
    169:     };
    170:     template<typename T>
    171:     struct mergePair: public std::binary_function<std::pair<T,T>, std::pair<T,T>, bool> {
    172:       std::pair<T,std::pair<T,T> > operator()(const std::pair<T,T>& x, const std::pair<T,T>& y) const {
    173:         return std::pair<T,std::pair<T,T> >(x.first, std::pair<T,T>(x.second, y.second));
    174:       }
    175:     };
    176:     template<typename _InputIterator1, typename _InputIterator2, typename _OutputIterator, typename _Compare, typename _Merge>
    177:     static _OutputIterator set_intersection_merge(_InputIterator1 __first1, _InputIterator1 __last1,
    178:                                            _InputIterator2 __first2, _InputIterator2 __last2,
    179:                                            _OutputIterator __result, _Compare __comp, _Merge __merge)
    180:     {
    181:       while(__first1 != __last1 && __first2 != __last2) {
    182:         if (__comp(*__first1, *__first2))
    183:           ++__first1;
    184:         else if (__comp(*__first2, *__first1))
    185:           ++__first2;
    186:         else
    187:         {
    188:           *__result = __merge(*__first1, *__first2);
    189:           ++__first1;
    190:           ++__first2;
    191:           ++__result;
    192:         }
    193:       }
    194:       return __result;
    195:     }
    196:   public:
    197:     template<typename Sequence, typename Renumbering, typename SendOverlap, typename RecvOverlap>
    198:     static void constructOverlap(const Sequence& points, Renumbering& renumbering, const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap) {
    199:       typedef typename SendOverlap::source_type point_type;
    200:       typedef std::pair<point_type,point_type>  pointPair;
    201:       typedef std::pair<point_type,pointPair>   pointTriple;
    202:       alloc_type allocator;
    203:       typename alloc_type::template rebind<point_type>::other point_allocator;
    204:       typename alloc_type::template rebind<pointPair>::other  pointPair_allocator;
    205:       const MPI_Comm comm     = sendOverlap->comm();
    206:       const int      commSize = sendOverlap->commSize();
    207:       const int      commRank = sendOverlap->commRank();
    208:       point_type    *sendBuf  = point_allocator.allocate(points.size()*2);
    209:       for(size_t i = 0; i < points.size()*2; ++i) {point_allocator.construct(sendBuf+i, point_type());}
    210:       int            size     = 0;
    211:       const int      debug    = sendOverlap->debug();
    212:       for(typename Sequence::const_iterator l_iter = points.begin(); l_iter != points.end(); ++l_iter) {
    213:         if (debug) {std::cout << "["<<commRank<<"]Send point["<<size<<"]: " << *l_iter << " " << renumbering[*l_iter] << std::endl;}
    214:         sendBuf[size++] = *l_iter;
    215:         sendBuf[size++] = renumbering[*l_iter];
    216:       }
    217:       int *sizes = allocator.allocate(commSize*3+2); // [size]   The number of points coming from each process
    218:       for(int i = 0; i < commSize*3+2; ++i) {allocator.construct(sizes+i, 0);}
    219:       int *offsets = sizes+commSize;                 // [size+1] Prefix sums for sizes
    220:       int *oldOffs = offsets+commSize+1;             // [size+1] Temporary storage
    221:       pointPair  *remotePoints = NULL;               // The points from each process
    222:       int        *remoteRanks  = NULL;               // The rank and number of overlap points of each process that overlaps another
    223:       int         numRemotePoints = 0;
    224:       int         numRemoteRanks  = 0;

    226:       // Change to Allgather() for the correct binning algorithm
    227:       MPI_Gather(&size, 1, MPI_INT, sizes, 1, MPI_INT, 0, comm);
    228:       if (commRank == 0) {
    229:         offsets[0] = 0;
    230:         for(int p = 1; p <= commSize; p++) {
    231:           offsets[p] = offsets[p-1] + sizes[p-1];
    232:         }
    233:         numRemotePoints = offsets[commSize];
    234:         remotePoints    = pointPair_allocator.allocate(numRemotePoints/2);
    235:         for(int i = 0; i < numRemotePoints/2; ++i) {pointPair_allocator.construct(remotePoints+i, pointPair());}
    236:       }
    237:       MPI_Gatherv(sendBuf, size, MPI_INT, remotePoints, sizes, offsets, MPI_INT, 0, comm);
    238:       for(size_t i = 0; i < points.size(); ++i) {point_allocator.destroy(sendBuf+i);}
    239:       point_allocator.deallocate(sendBuf, points.size());
    240:       std::map<int, std::map<int, std::set<pointTriple> > > overlapInfo; // Maps (p,q) to their set of overlap points

    242:       if (commRank == 0) {
    243:         for(int p = 0; p <= commSize; p++) {
    244:           offsets[p] /= 2;
    245:         }
    246:         for(int p = 0; p < commSize; p++) {
    247:           std::sort(&remotePoints[offsets[p]], &remotePoints[offsets[p+1]], lessPair<point_type>());
    248:         }
    249:         for(int p = 0; p <= commSize; p++) {
    250:           oldOffs[p] = offsets[p];
    251:         }
    252:         for(int p = 0; p < commSize; p++) {
    253:           for(int q = 0; q < commSize; q++) {
    254:             if (p == q) continue;
    255:             set_intersection_merge(&remotePoints[oldOffs[p]], &remotePoints[oldOffs[p+1]],
    256:                                    &remotePoints[oldOffs[q]], &remotePoints[oldOffs[q+1]],
    257:                                    std::insert_iterator<std::set<pointTriple> >(overlapInfo[p][q], overlapInfo[p][q].begin()),
    258:                                    lessPair<point_type>(), mergePair<point_type>());
    259:           }
    260:           sizes[p]     = overlapInfo[p].size()*2;
    261:           offsets[p+1] = offsets[p] + sizes[p];
    262:         }
    263:         numRemoteRanks = offsets[commSize];
    264:         remoteRanks    = allocator.allocate(numRemoteRanks);
    265:         for(int i = 0; i < numRemoteRanks; ++i) {allocator.construct(remoteRanks+i, 0);}
    266:         int     k = 0;
    267:         for(int p = 0; p < commSize; p++) {
    268:           for(typename std::map<int, std::set<pointTriple> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
    269:             remoteRanks[k*2]   = r_iter->first;
    270:             remoteRanks[k*2+1] = r_iter->second.size();
    271:             k++;
    272:           }
    273:         }
    274:       }
    275:       int numOverlaps;                          // The number of processes overlapping this process
    276:       MPI_Scatter(sizes, 1, MPI_INT, &numOverlaps, 1, MPI_INT, 0, comm);
    277:       int *overlapRanks = allocator.allocate(numOverlaps); // The rank and overlap size for each overlapping process
    278:       for(int i = 0; i < numOverlaps; ++i) {allocator.construct(overlapRanks+i, 0);}
    279:       MPI_Scatterv(remoteRanks, sizes, offsets, MPI_INT, overlapRanks, numOverlaps, MPI_INT, 0, comm);
    280:       point_type *sendPoints    = NULL;         // The points to send to each process
    281:       int         numSendPoints = 0;
    282:       if (commRank == 0) {
    283:         for(int p = 0, k = 0; p < commSize; p++) {
    284:           sizes[p] = 0;
    285:           for(int r = 0; r < (int) overlapInfo[p].size(); r++) {
    286:             sizes[p] += remoteRanks[k*2+1]*2;
    287:             k++;
    288:           }
    289:           offsets[p+1] = offsets[p] + sizes[p];
    290:         }
    291:         numSendPoints = offsets[commSize];
    292:         sendPoints    = point_allocator.allocate(numSendPoints*2);
    293:         for(int i = 0; i < numSendPoints*2; ++i) {point_allocator.construct(sendPoints+i, point_type());}
    294:         for(int p = 0, k = 0; p < commSize; p++) {
    295:           for(typename std::map<int, std::set<pointTriple> >::const_iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
    296:             int rank = r_iter->first;
    297:             for(typename std::set<pointTriple>::const_iterator p_iter = (overlapInfo[p][rank]).begin(); p_iter != (overlapInfo[p][rank]).end(); ++p_iter) {
    298:               sendPoints[k++] = p_iter->first;
    299:               sendPoints[k++] = p_iter->second.second;
    300:               if (debug) {std::cout << "["<<commRank<<"]Sending points " << p_iter->first << " " << p_iter->second.second << " to rank " << rank << std::endl;}
    301:             }
    302:           }
    303:         }
    304:       }
    305:       int numOverlapPoints = 0;
    306:       for(int r = 0; r < numOverlaps/2; r++) {
    307:         numOverlapPoints += overlapRanks[r*2+1];
    308:       }
    309:       point_type *overlapPoints = point_allocator.allocate(numOverlapPoints*2);
    310:       for(int i = 0; i < numOverlapPoints*2; ++i) {point_allocator.construct(overlapPoints+i, point_type());}
    311:       MPI_Scatterv(sendPoints, sizes, offsets, MPI_INT, overlapPoints, numOverlapPoints*2, MPI_INT, 0, comm);

    313:       for(int r = 0, k = 0; r < numOverlaps/2; r++) {
    314:         int rank = overlapRanks[r*2];

    316:         for(int p = 0; p < overlapRanks[r*2+1]; p++) {
    317:           point_type point       = overlapPoints[k++];
    318:           point_type remotePoint = overlapPoints[k++];

    320:           if (debug) {std::cout << "["<<commRank<<"]Matched up remote point " << remotePoint << "("<<point<<") to local " << renumbering[point] << std::endl;}
    321:           sendOverlap->addArrow(renumbering[point], rank, remotePoint);
    322:           recvOverlap->addArrow(rank, renumbering[point], remotePoint);
    323:         }
    324:       }

    326:       for(int i = 0; i < numOverlapPoints; ++i) {point_allocator.destroy(overlapPoints+i);}
    327:       point_allocator.deallocate(overlapPoints, numOverlapPoints);
    328:       for(int i = 0; i < numOverlaps; ++i) {allocator.destroy(overlapRanks+i);}
    329:       allocator.deallocate(overlapRanks, numOverlaps);
    330:       for(int i = 0; i < commSize*3+2; ++i) {allocator.destroy(sizes+i);}
    331:       allocator.deallocate(sizes, commSize*3+2);
    332:       if (commRank == 0) {
    333:         for(int i = 0; i < numRemoteRanks; ++i) {allocator.destroy(remoteRanks+i);}
    334:         allocator.deallocate(remoteRanks, numRemoteRanks);
    335:         for(int i = 0; i < numRemotePoints; ++i) {pointPair_allocator.destroy(remotePoints+i);}
    336:         pointPair_allocator.deallocate(remotePoints, numRemotePoints);
    337:         for(int i = 0; i < numSendPoints; ++i) {point_allocator.destroy(sendPoints+i);}
    338:         point_allocator.deallocate(sendPoints, numSendPoints);
    339:       }
    340:       // TODO: Rewrite above to use optimized construction
    341:       sendOverlap->assemble();
    342:       recvOverlap->assemble();
    343:       sendOverlap->assemblePoints();
    344:       recvOverlap->assemblePoints();
    345:     }
    346:   };
    347:   namespace Pullback {
    348:     class SimpleCopy {
    349:     public:
    350:       // Copy the overlap section to the related processes
    351:       //   This version is for Constant sections, meaning the same, single value over all points
    352:       template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
    353:       static void copyConstant(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
    354:         MPIMover<char>                             pMover(sendSection->comm(), sendSection->debug());
    355:         MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), sendSection->debug());
    356:         std::map<int, char *>                      sendPoints;
    357:         std::map<int, char *>                      recvPoints;
    358:         typename SendSection::alloc_type::template rebind<char>::other sendAllocator;
    359:         typename RecvSection::alloc_type::template rebind<char>::other recvAllocator;

    361:         const typename SendOverlap::baseSequence::iterator sBegin  = sendOverlap->baseBegin();
    362:         const typename SendOverlap::baseSequence::iterator sEnd    = sendOverlap->baseEnd();
    363:         const typename SendSection::value_type            *sValues = sendSection->restrictSpace();

    365:         for(typename SendOverlap::baseSequence::iterator r_iter = sBegin; r_iter != sEnd; ++r_iter) {
    366:           const int                                          pSize  = sendOverlap->getConeSize(*r_iter);
    367:           const typename SendOverlap::coneSequence::iterator pBegin = sendOverlap->coneBegin(*r_iter);
    368:           const typename SendOverlap::coneSequence::iterator pEnd   = sendOverlap->coneEnd(*r_iter);
    369:           char                                              *v      = sendAllocator.allocate(pSize);
    370:           int                                                k      = 0;

    372:           for(int i = 0; i < pSize; ++i) {sendAllocator.construct(v+i, 0);}
    373:           for(typename SendOverlap::coneSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter, ++k) {
    374:             v[k] = (char) sendSection->hasPoint(*p_iter);
    375:           }
    376:           sendPoints[*r_iter] = v;
    377:           pMover.send(*r_iter, pSize, sendPoints[*r_iter]);
    378:           vMover.send(*r_iter, 2, sValues);
    379:         }
    380:         const typename RecvOverlap::capSequence::iterator rBegin  = recvOverlap->capBegin();
    381:         const typename RecvOverlap::capSequence::iterator rEnd    = recvOverlap->capEnd();
    382:         const typename RecvSection::value_type           *rValues = recvSection->restrictSpace();

    384:         for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
    385:           const int pSize = recvOverlap->getSupportSize(*r_iter);
    386:           char     *v     = recvAllocator.allocate(pSize);

    388:           for(int i = 0; i < pSize; ++i) {recvAllocator.construct(v+i, 0);}
    389:           recvPoints[*r_iter] = v;
    390:           pMover.recv(*r_iter, pSize, recvPoints[*r_iter]);
    391:           vMover.recv(*r_iter, 2, rValues);
    392:         }
    393:         pMover.start();
    394:         pMover.end();
    395:         vMover.start();
    396:         vMover.end();
    397:         for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
    398:           const typename RecvOverlap::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
    399:           const typename RecvOverlap::supportSequence::iterator pEnd   = recvOverlap->supportEnd(*r_iter);
    400:           const char                                           *v      = recvPoints[*r_iter];
    401:           int                                                   k      = 0;

    403:           for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter, ++k) {
    404:             if (v[k]) {recvSection->addPoint(typename RecvSection::point_type(*r_iter, s_iter.color()));}
    405:           }
    406:         }
    407:         for(typename SendOverlap::baseSequence::iterator r_iter = sBegin; r_iter != sEnd; ++r_iter) {
    408:           sendAllocator.deallocate(sendPoints[*r_iter], sendOverlap->getConeSize(*r_iter));
    409:         }
    410:         for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
    411:           recvAllocator.deallocate(recvPoints[*r_iter], recvOverlap->getSupportSize(*r_iter));
    412:         }
    413:       }
    414:       // Specialize to ArrowSection
    415:       template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
    416:       static void copyConstantArrow(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
    417:         MPIMover<char>                             pMover(sendSection->comm(), sendSection->debug());
    418:         MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), sendSection->debug());
    419:         std::map<int, char *>                      sendPoints;
    420:         std::map<int, char *>                      recvPoints;
    421:         typename SendSection::alloc_type::template rebind<char>::other sendAllocator;
    422:         typename RecvSection::alloc_type::template rebind<char>::other recvAllocator;

    424:         const Obj<typename SendOverlap::traits::baseSequence>      sRanks  = sendOverlap->base();
    425:         const typename SendOverlap::traits::baseSequence::iterator sEnd    = sRanks->end();
    426:         const typename SendSection::value_type                    *sValues = sendSection->restrictSpace();

    428:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
    429:           const Obj<typename SendOverlap::coneSequence>&     points = sendOverlap->cone(*r_iter);
    430:           const int                                          pSize  = sendOverlap->getConeSize(*r_iter);
    431:           const typename SendOverlap::coneSequence::iterator pBegin = points->begin();
    432:           const typename SendOverlap::coneSequence::iterator pEnd   = points->end();
    433:           char                                              *v      = sendAllocator.allocate(pSize*pSize);
    434:           int                                                k      = 0;

    436:           for(size_t i = 0; i < pSize*pSize; ++i) {sendAllocator.construct(v+i, 0);}
    437:           for(typename SendOverlap::coneSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
    438:             for(typename SendOverlap::coneSequence::iterator q_iter = pBegin; q_iter != pEnd; ++q_iter, ++k) {
    439:               v[k] = (char) sendSection->hasPoint(typename SendSection::point_type(*p_iter,*q_iter));
    440:             }
    441:           }
    442:           sendPoints[*r_iter] = v;
    443:           pMover.send(*r_iter, pSize*pSize, sendPoints[*r_iter]);
    444:           vMover.send(*r_iter, 2, sValues);
    445:         }
    446:         const Obj<typename RecvOverlap::traits::capSequence>      rRanks  = recvOverlap->cap();
    447:         const typename RecvOverlap::traits::capSequence::iterator rEnd    = rRanks->end();
    448:         const typename RecvSection::value_type                   *rValues = recvSection->restrictSpace();

    450:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
    451:           const Obj<typename RecvOverlap::traits::supportSequence>& points = recvOverlap->support(*r_iter);
    452:           const int                                                 pSize  = recvOverlap->getSupportSize(*r_iter);
    453:           char                                                     *v      = recvAllocator.allocate(pSize*pSize);

    455:           for(size_t i = 0; i < pSize*pSize; ++i) {recvAllocator.construct(v+i, 0);}
    456:           recvPoints[*r_iter] = v;
    457:           pMover.recv(*r_iter, pSize*pSize, recvPoints[*r_iter]);
    458:           vMover.recv(*r_iter, 2, rValues);
    459:         }
    460:         pMover.start();
    461:         pMover.end();
    462:         vMover.start();
    463:         vMover.end();
    464:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
    465:           const Obj<typename RecvOverlap::traits::supportSequence>&     points = recvOverlap->support(*r_iter);
    466:           const typename RecvOverlap::traits::supportSequence::iterator pBegin = points->begin();
    467:           const typename RecvOverlap::traits::supportSequence::iterator pEnd   = points->end();
    468:           const char                                                   *v      = recvPoints[*r_iter];
    469:           int                                                           k      = 0;

    471:           for(typename RecvOverlap::traits::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
    472:             for(typename RecvOverlap::traits::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter, ++k) {
    473:               if (v[k]) {recvSection->addPoint(typename RecvSection::point_type(s_iter.color(),t_iter.color()));}
    474:             }
    475:           }
    476:         }
    477:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
    478:           sendAllocator.deallocate(sendPoints[*r_iter], sendOverlap->getConeSize(*r_iter)*sendOverlap->getConeSize(*r_iter));
    479:         }
    480:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
    481:           recvAllocator.deallocate(recvPoints[*r_iter], recvOverlap->getSupportSize(*r_iter)*recvOverlap->getSupportSize(*r_iter));
    482:         }
    483:       }
    484:       // Copy the overlap section to the related processes
    485:       //   This version is for IConstant sections, meaning the same, single value over all points
    486:       template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
    487:       static void copyIConstant(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
    488:         MPIMover<typename SendSection::point_type> pMover(sendSection->comm(), sendSection->debug());
    489:         MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), sendSection->debug());
    490:         std::map<int, typename SendSection::point_type *> sendPoints;
    491:         std::map<int, typename SendSection::point_type *> recvPoints;
    492:         typename SendSection::alloc_type::template rebind<typename SendSection::point_type>::other sendAllocator;
    493:         typename RecvSection::alloc_type::template rebind<typename SendSection::point_type>::other recvAllocator;

    495:         const Obj<typename SendOverlap::baseSequence>      sRanks  = sendOverlap->base();
    496:         const typename SendOverlap::baseSequence::iterator sEnd    = sRanks->end();
    497:         const typename SendSection::value_type            *sValues = sendSection->restrictSpace();

    499:         for(typename SendOverlap::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
    500:           typename SendSection::point_type *v = sendAllocator.allocate(2);

    502:           for(size_t i = 0; i < 2; ++i) {sendAllocator.construct(v+i, 0);}
    503:           v[0] = sendSection->getChart().min();
    504:           v[1] = sendSection->getChart().max();
    505:           sendPoints[*r_iter] = v;
    506:           pMover.send(*r_iter, 2, sendPoints[*r_iter]);
    507:           vMover.send(*r_iter, 2, sValues);
    508:           ///std::cout << "["<<sendOverlap->commRank()<<"]Sending chart (" << v[0] << ", " << v[1] << ") with values (" << sValues[0] << ", " << sValues[1] << ") to process " << *r_iter << std::endl;
    509:         }
    510:         const Obj<typename RecvOverlap::capSequence>      rRanks  = recvOverlap->cap();
    511:         const typename RecvOverlap::capSequence::iterator rEnd    = rRanks->end();
    512:         const typename RecvSection::value_type           *rValues = recvSection->restrictSpace();

    514:         for(typename RecvOverlap::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
    515:           typename SendSection::point_type *v = recvAllocator.allocate(2);

    517:           for(size_t i = 0; i < 2; ++i) {recvAllocator.construct(v+i, 0);}
    518:           recvPoints[*r_iter] = v;
    519:           pMover.recv(*r_iter, 2, recvPoints[*r_iter]);
    520:           vMover.recv(*r_iter, 2, rValues);
    521:         }
    522:         pMover.start();
    523:         pMover.end();
    524:         vMover.start();
    525:         vMover.end();

    527:         typename SendSection::point_type min = -1;
    528:         typename SendSection::point_type max = -1;

    530:         for(typename RecvOverlap::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
    531:           const typename RecvSection::point_type *v = recvPoints[*r_iter];
    532:           typename SendSection::point_type        newMin = v[0];
    533:           typename SendSection::point_type        newMax = v[1]-1;
    534:           //int                                     pSize  = 0;

    536:           ///std::cout << "["<<recvOverlap->commRank()<<"]Received chart (" << v[0] << ", " << v[1] << ") from process " << *r_iter << std::endl;
    537: #if 0
    538:           // Translate to local numbering
    539:           if (recvOverlap->support(*r_iter)->size()) {
    540:             while(!pSize) {
    541:               const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter, newMin);
    542:               pSize = points->size();
    543:               if (pSize) {
    544:                 newMin = *points->begin();
    545:               } else {
    546:                 newMin++;
    547:               }
    548:             }
    549:             pSize  = 0;
    550:             while(!pSize) {
    551:               const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter, newMax);
    552:               pSize = points->size();
    553:               if (pSize) {
    554:                 newMax = *points->begin();
    555:               } else {
    556:                 newMax--;
    557:               }
    558:             }
    559:           }
    560:           std::cout << "["<<recvOverlap->commRank()<<"]Translated to chart (" << newMin << ", " << newMax+1 << ") from process " << *r_iter << std::endl;
    561: #endif
    562:           // Update chart
    563:           if (min < 0) {
    564:             min = newMin;
    565:             max = newMax+1;
    566:           } else {
    567:             min = std::min(min, newMin);
    568:             max = std::max(max, (typename SendSection::point_type) (newMax+1));
    569:           }
    570:         }
    571:         if (!rRanks->size()) {min = max = 0;}
    572:         recvSection->setChart(typename RecvSection::chart_type(min, max));
    573:         for(typename SendOverlap::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
    574:           sendAllocator.deallocate(sendPoints[*r_iter], 2);
    575:         }
    576:         for(typename RecvOverlap::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
    577:           recvAllocator.deallocate(recvPoints[*r_iter], 2);
    578:         }
    579:       }
    580:       // Copy the overlap section to the related processes
    581:       //   This version is for different sections, possibly with different data types
    582:       // TODO: Can cache MPIMover objects (like a VecScatter)
    583:       template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
    584:       static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection, const MPI_Datatype datatype = MPI_DATATYPE_NULL) {
    585:         typedef std::pair<int, typename SendSection::value_type *> allocPair;
    586:         typedef typename RecvSection::point_type                   recv_point_type;
    587:         const Obj<typename SendSection::atlas_type>& sendAtlas = sendSection->getAtlas();
    588:         const Obj<typename RecvSection::atlas_type>& recvAtlas = recvSection->getAtlas();
    589:         MPIMover<typename SendSection::value_type>   vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
    590:         std::map<int, allocPair>                     sendValues;
    591:         std::map<int, allocPair>                     recvValues;
    592:         typename SendSection::alloc_type             sendAllocator;
    593:         typename RecvSection::alloc_type             recvAllocator;

    595:         copy(sendOverlap, recvOverlap, sendAtlas, recvAtlas);
    596:         const typename SendOverlap::baseSequence::iterator sBegin = sendOverlap->baseBegin();
    597:         const typename SendOverlap::baseSequence::iterator sEnd   = sendOverlap->baseEnd();

    599:         // TODO: This should be const_iterator, but Sifter sucks
    600:         for(typename SendOverlap::baseSequence::iterator r_iter = sBegin; r_iter != sEnd; ++r_iter) {
    601:           const typename SendOverlap::coneSequence::iterator pBegin    = sendOverlap->coneBegin(*r_iter);
    602:           const typename SendOverlap::coneSequence::iterator pEnd      = sendOverlap->coneEnd(*r_iter);
    603:           const int                                          numPoints = sendOverlap->getConeSize(*r_iter);
    604:           std::valarray<typename SendOverlap::source_type>   sortedPoints(numPoints);
    605:           int                                                numVals   = 0, p = 0;

    607:           // TODO: This should be const_iterator, but Sifter sucks
    608:           for(typename SendOverlap::coneSequence::iterator c_iter = pBegin; c_iter != pEnd; ++c_iter, ++p) {
    609:             numVals += sendSection->getFiberDimension(*c_iter);
    610:             sortedPoints[p] = *c_iter;
    611:           }
    612:           typename SendSection::value_type *v = sendAllocator.allocate(numVals);
    613:           int                               k = 0;

    615:           std::sort(&sortedPoints[0], &sortedPoints[numPoints]);
    616:           for(int i = 0; i < numVals; ++i) {sendAllocator.construct(v+i, 0);}
    617:           for(p = 0; p < numPoints; ++p) {
    618:             const typename SendSection::value_type *vals = sendSection->restrictPoint(sortedPoints[p]);

    620:             for(int i = 0; i < sendSection->getFiberDimension(sortedPoints[p]); ++i, ++k) v[k] = vals[i];
    621:           }
    622:           sendValues[*r_iter] = allocPair(numVals, v);
    623:           vMover.send(*r_iter, numVals, v);
    624:         }
    625:         const typename RecvOverlap::capSequence::iterator rBegin = recvOverlap->capBegin();
    626:         const typename RecvOverlap::capSequence::iterator rEnd   = recvOverlap->capEnd();

    628:         recvSection->allocatePoint();
    629:         // TODO: This should be const_iterator, but Sifter sucks
    630:         int maxVals = 0;
    631:         for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
    632:           const typename RecvOverlap::supportSequence::iterator pBegin  = recvOverlap->supportBegin(*r_iter);
    633:           const typename RecvOverlap::supportSequence::iterator pEnd    = recvOverlap->supportEnd(*r_iter);
    634:           int                                                   numVals = 0;

    636:           // TODO: This should be const_iterator, but Sifter sucks
    637:           for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
    638:             numVals += recvSection->getFiberDimension(recv_point_type(*r_iter, s_iter.color()));
    639:           }
    640:           typename SendSection::value_type *v = sendAllocator.allocate(numVals);

    642:           for(int i = 0; i < numVals; ++i) {sendAllocator.construct(v+i, 0);}
    643:           recvValues[*r_iter] = allocPair(numVals, v);
    644:           vMover.recv(*r_iter, numVals, v);
    645:           maxVals = std::max(maxVals, numVals);
    646:         }
    647:         vMover.start();
    648:         vMover.end();
    649:         typename RecvSection::value_type *convertedValues = recvAllocator.allocate(maxVals);
    650:         for(int i = 0; i < maxVals; ++i) {recvAllocator.construct(convertedValues+i, 0);}
    651:         for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
    652:           const typename RecvOverlap::supportSequence::iterator pBegin    = recvOverlap->supportBegin(*r_iter);
    653:           const typename RecvOverlap::supportSequence::iterator pEnd      = recvOverlap->supportEnd(*r_iter);
    654:           const int                                             numPoints = recvOverlap->getSupportSize(*r_iter);
    655:           std::valarray<typename RecvOverlap::color_type>       sortedPoints(numPoints);
    656:           typename SendSection::value_type                     *v       = recvValues[*r_iter].second;
    657:           const int                                             numVals = recvValues[*r_iter].first;
    658:           int                                                   k       = 0, p = 0;

    660:           for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter, ++p) {
    661:             sortedPoints[p] = s_iter.color();
    662:           }
    663:           std::sort(&sortedPoints[0], &sortedPoints[numPoints]);

    665:           //for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter) {
    666:           for(p = 0; p < numPoints; ++p) {
    667:             const int size = recvSection->getFiberDimension(recv_point_type(*r_iter, sortedPoints[p]));

    669:             for(int i = 0; i < size; ++i) {convertedValues[i] = (typename RecvSection::value_type) v[k+i];}
    670:             if (size) {recvSection->updatePoint(recv_point_type(*r_iter, sortedPoints[p]), convertedValues);}
    671:             k += size;
    672:           }
    673:           assert(k == numVals);
    674:           for(int i = 0; i < numVals; ++i) {sendAllocator.destroy(v+i);}
    675:         }
    676:         for(int i = 0; i < maxVals; ++i) {recvAllocator.destroy(convertedValues+i);}
    677:         recvAllocator.deallocate(convertedValues, maxVals);
    678:         for(typename SendOverlap::baseSequence::iterator r_iter = sBegin; r_iter != sEnd; ++r_iter) {
    679:           typename SendSection::value_type *v       = sendValues[*r_iter].second;
    680:           const int                         numVals = sendValues[*r_iter].first;

    682:           for(int i = 0; i < numVals; ++i) {sendAllocator.destroy(v+i);}
    683:           sendAllocator.deallocate(v, numVals);
    684:         }
    685:         for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
    686:           typename SendSection::value_type *v       = recvValues[*r_iter].second;
    687:           const int                         numVals = recvValues[*r_iter].first;

    689:           for(int i = 0; i < numVals; ++i) {sendAllocator.destroy(v+i);}
    690:           sendAllocator.deallocate(v, numVals);
    691:         }
    692:         //recvSection->view("After copy");
    693:       }
    694:       // Copy the overlap section to the related processes
    695:       //   This version is for sections with the same type
    696:       template<typename SendOverlap, typename RecvOverlap, typename Section>
    697:       static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& sendSection, const Obj<Section>& recvSection, const MPI_Datatype datatype = MPI_DATATYPE_NULL) {
    698:         typedef std::pair<int, typename Section::value_type *> allocPair;
    699:         const Obj<typename Section::atlas_type>& sendAtlas = sendSection->getAtlas();
    700:         const Obj<typename Section::atlas_type>& recvAtlas = recvSection->getAtlas();
    701:         MPIMover<typename Section::value_type>   vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
    702:         std::map<int, allocPair>                 sendValues;
    703:         std::map<int, allocPair>                 recvValues;
    704:         typename Section::alloc_type             allocator;

    706:         ///sendAtlas->view("Send Atlas in same type copy()");
    707:         copy(sendOverlap, recvOverlap, sendAtlas, recvAtlas);
    708:         ///recvAtlas->view("Recv Atlas after same type copy()");
    709:         const Obj<typename SendOverlap::traits::baseSequence>      sRanks = sendOverlap->base();
    710:         const typename SendOverlap::traits::baseSequence::iterator sEnd   = sRanks->end();

    712:         // TODO: This should be const_iterator, but Sifter sucks
    713:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
    714:           const Obj<typename SendOverlap::coneSequence>&     points    = sendOverlap->cone(*r_iter);
    715:           const typename SendOverlap::coneSequence::iterator pEnd      = points->end();
    716:           const int                                          numPoints = sendOverlap->getConeSize(*r_iter);
    717:           std::valarray<typename SendOverlap::source_type>   sortedPoints(numPoints);
    718:           int                                                numVals   = 0, p = 0;

    720:           // TODO: This should be const_iterator, but Sifter sucks
    721:           for(typename SendOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != pEnd; ++c_iter, ++p) {
    722:             numVals += sendSection->getFiberDimension(*c_iter);
    723:             sortedPoints[p] = *c_iter;
    724:           }
    725:           typename Section::value_type *v = allocator.allocate(numVals);
    726:           int                           k = 0;

    728:           std::sort(&sortedPoints[0], &sortedPoints[numPoints]);
    729:           for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
    730:           //for(typename SendOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != pEnd; ++c_iter) {
    731:           for(p = 0; p < numPoints; ++p) {
    732:             const typename Section::value_type *vals = sendSection->restrictPoint(sortedPoints[p]);
    733:             const int                           dim  = sendSection->getFiberDimension(sortedPoints[p]);

    735:             for(int i = 0; i < dim; ++i, ++k) v[k] = vals[i];
    736:           }
    737:           sendValues[*r_iter] = allocPair(numVals, v);
    738:           vMover.send(*r_iter, numVals, v);
    739:         }
    740:         const Obj<typename RecvOverlap::traits::capSequence>      rRanks = recvOverlap->cap();
    741:         const typename RecvOverlap::traits::capSequence::iterator rEnd   = rRanks->end();

    743:         recvSection->allocatePoint();
    744:         ///recvSection->view("Recv Section after same type copy() allocation");
    745:         // TODO: This should be const_iterator, but Sifter sucks
    746:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
    747:           const Obj<typename RecvOverlap::supportSequence>&     points  = recvOverlap->support(*r_iter);
    748:           const typename RecvOverlap::supportSequence::iterator pEnd    = points->end();
    749:           int                                                   numVals = 0;

    751:           // TODO: This should be const_iterator, but Sifter sucks
    752:           for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter) {
    753:             numVals += recvSection->getFiberDimension(s_iter.color());
    754:           }
    755:           typename Section::value_type *v = allocator.allocate(numVals);

    757:           recvValues[*r_iter] = allocPair(numVals, v);
    758:           for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
    759:           vMover.recv(*r_iter, numVals, v);
    760:         }
    761:         vMover.start();
    762:         vMover.end();
    763:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
    764:           const Obj<typename RecvOverlap::supportSequence>&     points    = recvOverlap->support(*r_iter);
    765:           const typename RecvOverlap::supportSequence::iterator pEnd      = points->end();
    766:           const int                                             numPoints = recvOverlap->getSupportSize(*r_iter);
    767:           std::valarray<typename RecvOverlap::color_type>       sortedPoints(numPoints);
    768:           typename Section::value_type                         *v         = recvValues[*r_iter].second;
    769:           const int                                             numVals   = recvValues[*r_iter].first;
    770:           int                                                   k         = 0, p = 0;

    772:           for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter, ++p) {
    773:             sortedPoints[p] = s_iter.color();
    774:           }
    775:           std::sort(&sortedPoints[0], &sortedPoints[numPoints]);

    777:           //for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter) {
    778:           for(p = 0; p < numPoints; ++p) {
    779:             const int size = recvSection->getFiberDimension(sortedPoints[p]);

    781:             if (size) {recvSection->updatePoint(sortedPoints[p], &v[k]);}
    782:             k += size;
    783:           }
    784:           assert(k == numVals);
    785:           for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
    786:           allocator.deallocate(v, numVals);
    787:         }
    788:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
    789:           typename Section::value_type *v       = sendValues[*r_iter].second;
    790:           const int                     numVals = sendValues[*r_iter].first;

    792:           for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
    793:           allocator.deallocate(v, numVals);
    794:         }
    795:         //recvSection->view("After copy");
    796:       }
    797:       // Specialize to ArrowSection
    798:       template<typename SendOverlap, typename RecvOverlap>
    799:       static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<UniformSection<MinimalArrow<int,int>,int> >& sendSection, const Obj<UniformSection<MinimalArrow<int,int>,int> >& recvSection, const MPI_Datatype datatype = MPI_DATATYPE_NULL) {
    800:         typedef UniformSection<MinimalArrow<int,int>,int>      Section;
    801:         typedef std::pair<int, typename Section::value_type *> allocPair;
    802:         const Obj<typename Section::atlas_type>& sendAtlas = sendSection->getAtlas();
    803:         const Obj<typename Section::atlas_type>& recvAtlas = recvSection->getAtlas();
    804:         MPIMover<typename Section::value_type>   vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
    805:         std::map<int, allocPair>                 sendValues;
    806:         std::map<int, allocPair>                 recvValues;
    807:         typename Section::alloc_type             allocator;

    809:         copy(sendOverlap, recvOverlap, sendAtlas, recvAtlas);
    810:         const Obj<typename SendOverlap::traits::baseSequence>      sRanks = sendOverlap->base();
    811:         const typename SendOverlap::traits::baseSequence::iterator sEnd   = sRanks->end();

    813:         // TODO: This should be const_iterator, but Sifter sucks
    814:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
    815:           const Obj<typename SendOverlap::coneSequence>&     points  = sendOverlap->cone(*r_iter);
    816:           const typename SendOverlap::coneSequence::iterator pBegin  = points->begin();
    817:           const typename SendOverlap::coneSequence::iterator pEnd    = points->end();
    818:           int                                                numVals = 0;

    820:           // TODO: This should be const_iterator, but Sifter sucks
    821:           for(typename SendOverlap::coneSequence::iterator c_iter = pBegin; c_iter != pEnd; ++c_iter) {
    822:             for(typename SendOverlap::coneSequence::iterator d_iter = pBegin; d_iter != pEnd; ++d_iter) {
    823:               numVals += sendSection->getFiberDimension(typename Section::point_type(*c_iter, *d_iter));
    824:             }
    825:           }
    826:           typename Section::value_type *v = allocator.allocate(numVals);
    827:           int                           k = 0;

    829:           for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
    830:           for(typename SendOverlap::coneSequence::iterator c_iter = pBegin; c_iter != pEnd; ++c_iter) {
    831:             for(typename SendOverlap::coneSequence::iterator d_iter = pBegin; d_iter != pEnd; ++d_iter) {
    832:               const typename Section::point_type  arrow(*c_iter, *d_iter);
    833:               const typename Section::value_type *vals = sendSection->restrictPoint(arrow);
    834:               const int                           dim  = sendSection->getFiberDimension(arrow);

    836:               for(int i = 0; i < dim; ++i, ++k) v[k] = vals[i];
    837:             }
    838:           }
    839:           sendValues[*r_iter] = allocPair(numVals, v);
    840:           vMover.send(*r_iter, numVals, v);
    841:         }
    842:         const Obj<typename RecvOverlap::traits::capSequence>      rRanks = recvOverlap->cap();
    843:         const typename RecvOverlap::traits::capSequence::iterator rEnd   = rRanks->end();

    845:         recvSection->allocatePoint();
    846:         // TODO: This should be const_iterator, but Sifter sucks
    847:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
    848:           const Obj<typename RecvOverlap::supportSequence>&     points  = recvOverlap->support(*r_iter);
    849:           const typename RecvOverlap::supportSequence::iterator pBegin  = points->begin();
    850:           const typename RecvOverlap::supportSequence::iterator pEnd    = points->end();
    851:           int                                                   numVals = 0;

    853:           // TODO: This should be const_iterator, but Sifter sucks
    854:           for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
    855:             for(typename RecvOverlap::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter) {
    856:               numVals += recvSection->getFiberDimension(typename Section::point_type(s_iter.color(), t_iter.color()));
    857:             }
    858:           }
    859:           typename Section::value_type *v = allocator.allocate(numVals);

    861:           recvValues[*r_iter] = allocPair(numVals, v);
    862:           for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
    863:           vMover.recv(*r_iter, numVals, v);
    864:         }
    865:         vMover.start();
    866:         vMover.end();
    867:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
    868:           const Obj<typename RecvOverlap::supportSequence>&     points  = recvOverlap->support(*r_iter);
    869:           const typename RecvOverlap::supportSequence::iterator pBegin  = points->begin();
    870:           const typename RecvOverlap::supportSequence::iterator pEnd    = points->end();
    871:           typename Section::value_type                         *v       = recvValues[*r_iter].second;
    872:           const int                                             numVals = recvValues[*r_iter].first;
    873:           int                                                   k       = 0;

    875:           for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
    876:             for(typename RecvOverlap::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter) {
    877:               const typename Section::point_type arrow(s_iter.color(), t_iter.color());
    878:               const int size = recvSection->getFiberDimension(arrow);

    880:               if (size) {recvSection->updatePoint(arrow, &v[k]);}
    881:               k += size;
    882:             }
    883:           }
    884:           for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
    885:           allocator.deallocate(v, numVals);
    886:         }
    887:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
    888:           typename Section::value_type *v       = sendValues[*r_iter].second;
    889:           const int                     numVals = sendValues[*r_iter].first;

    891:           for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
    892:           allocator.deallocate(v, numVals);
    893:         }
    894:         //recvSection->view("After copy");
    895:       }
    896:       // Specialize to a ConstantSection
    897: #if 0
    898:       template<typename SendOverlap, typename RecvOverlap, typename Value>
    899:       static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<ConstantSection<typename SendOverlap::source_type, Value> >& sendSection, const Obj<ConstantSection<typename SendOverlap::source_type, Value> >& recvSection) {
    900:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
    901:       };
    902:       template<typename SendOverlap, typename RecvOverlap, typename Value>
    903:       static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<IConstantSection<typename SendOverlap::source_type, Value> >& sendSection, const Obj<ConstantSection<typename SendOverlap::source_type, Value> >& recvSection) {
    904:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
    905:       };
    906: #else
    907:       template<typename SendOverlap, typename RecvOverlap, typename Value>
    908:       static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<ConstantSection<typename SendOverlap::source_type, Value> >& sendSection, const Obj<ConstantSection<ALE::Pair<int, typename SendOverlap::source_type>, Value> >& recvSection) {
    909:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
    910:       }
    911:       template<typename SendOverlap, typename RecvOverlap, typename Value>
    912:       static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<IConstantSection<typename SendOverlap::source_type, Value> >& sendSection, const Obj<ConstantSection<ALE::Pair<int, typename SendOverlap::source_type>, Value> >& recvSection) {
    913:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
    914:       }
    915: #endif
    916:       // Specialize to an IConstantSection
    917:       template<typename SendOverlap, typename RecvOverlap, typename Value>
    918:       static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<IConstantSection<typename SendOverlap::source_type, Value> >& sendSection, const Obj<IConstantSection<typename SendOverlap::source_type, Value> >& recvSection) {
    919:         // Why doesn't this work?
    920:         //   This supposed to be a copy, BUT filtered through the sendOverlap
    921:         //   However, an IConstant section does not allow filtration of its
    922:         //   chart. Therefore, you end up with either
    923:         //
    924:         //   a) Too many items in the chart, copied from the remote sendSection
    925:         //   b) A chart mapped to the local numbering, which we do not want
    926:         copyIConstant(sendOverlap, recvOverlap, sendSection, recvSection);
    927:       }
    928:       // Specialize to an BaseSection/ConstantSection pair
    929: #if 0
    930:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
    931:       static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<BaseSection<Sieve_> >& sendSection, const Obj<ConstantSection<typename SendOverlap::source_type, int> >& recvSection) {
    932:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
    933:       };
    934: #else
    935:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
    936:       static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<BaseSection<Sieve_> >& sendSection, const Obj<ConstantSection<ALE::Pair<int, typename SendOverlap::source_type>, int> >& recvSection) {
    937:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
    938:       }
    939: #endif
    940:       // Specialize to an BaseSectionV/ConstantSection pair
    941: #if 0
    942:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
    943:       static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<BaseSectionV<Sieve_> >& sendSection, const Obj<ConstantSection<typename SendOverlap::source_type, int> >& recvSection) {
    944:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
    945:       };
    946: #else
    947:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
    948:       static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<BaseSectionV<Sieve_> >& sendSection, const Obj<ConstantSection<ALE::Pair<int, typename SendOverlap::source_type>, int> >& recvSection) {
    949:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
    950:       }
    951: #endif
    952:       // Specialize to an LabelBaseSection/ConstantSection pair
    953: #if 0
    954:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_, typename Label_>
    955:       static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<LabelBaseSection<Sieve_, Label_> >& sendSection, const Obj<ConstantSection<typename SendOverlap::source_type, int> >& recvSection) {
    956:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
    957:       };
    958: #else
    959:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_, typename Label_>
    960:       static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<LabelBaseSection<Sieve_, Label_> >& sendSection, const Obj<ConstantSection<ALE::Pair<int, typename SendOverlap::source_type>, int> >& recvSection) {
    961:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
    962:       }
    963: #endif
    964:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_, typename Label_>
    965:       static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<LabelBaseSectionV<Sieve_, Label_> >& sendSection, const Obj<ConstantSection<ALE::Pair<int, typename SendOverlap::source_type>, int> >& recvSection) {
    966:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
    967:       }
    968:       // Specialize to a ConstantSection for ArrowSection
    969:       template<typename SendOverlap, typename RecvOverlap, typename Value>
    970:       static void copy(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<ConstantSection<MinimalArrow<typename SendOverlap::source_type,typename SendOverlap::source_type>, Value> >& sendSection, const Obj<ConstantSection<MinimalArrow<typename SendOverlap::source_type,typename SendOverlap::source_type>, Value> >& recvSection) {
    971:         copyConstantArrow(sendOverlap, recvOverlap, sendSection, recvSection);
    972:       }
    973:     };
    974:     class BinaryFusion {
    975:     public:
    976:       template<typename OverlapSection, typename RecvOverlap, typename Section, typename Function>
    977:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section, Function binaryOp) {
    978:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
    979:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

    981:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
    982:           // TODO: This must become a more general iterator over colors
    983:           const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
    984:           // Just taking the first value
    985:           const typename Section::point_type&        localPoint    = *p_iter;
    986:           const typename OverlapSection::point_type& remotePoint   = points->begin().color();
    987:           const typename OverlapSection::value_type *overlapValues = overlapSection->restrictPoint(remotePoint);
    988:           const typename Section::value_type        *localValues   = section->restrictPoint(localPoint);
    989:           const int                                  dim           = section->getFiberDimension(localPoint);
    990:           // TODO: optimize allocation
    991:           typename Section::value_type              *values        = new typename Section::value_type[dim];

    993:           for(int d = 0; d < dim; ++d) {
    994:             values[d] = binaryOp(overlapValues[d], localValues[d]);
    995:           }
    996:           section->updatePoint(localPoint, values);
    997:           delete [] values;
    998:         }
    999:       }
    1000:     };
    1001:     class ReplacementBinaryFusion {
    1002:     public:
    1003:       template<typename OverlapSection, typename RecvOverlap, typename Section>
    1004:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
    1005:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
    1006:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

    1008:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
    1009:           // TODO: This must become a more general iterator over colors
    1010:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
    1011:           // Just taking the first value
    1012:           const typename Section::point_type&            localPoint  = *p_iter;
    1013:           const typename OverlapSection::point_type&     remotePoint = points->begin().color();

    1015:           section->update(localPoint, overlapSection->restrictPoint(remotePoint));
    1016:         }
    1017:       }
    1018:     };
    1019:     class AdditiveBinaryFusion {
    1020:     public:
    1021:       template<typename OverlapSection, typename RecvOverlap, typename Section>
    1022:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
    1023:         typedef typename Section::point_type        point_type;
    1024:         typedef typename OverlapSection::point_type overlap_point_type;
    1025:         const typename RecvOverlap::capSequence::iterator rBegin = recvOverlap->capBegin();
    1026:         const typename RecvOverlap::capSequence::iterator rEnd   = recvOverlap->capEnd();

    1028:         for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
    1029:           const int                                             rank    = *r_iter;
    1030:           const typename RecvOverlap::supportSequence::iterator pBegin  = recvOverlap->supportBegin(*r_iter);
    1031:           const typename RecvOverlap::supportSequence::iterator pEnd    = recvOverlap->supportEnd(*r_iter);

    1033:           for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
    1034:             const point_type& localPoint  = *p_iter;
    1035:             const point_type& remotePoint = p_iter.color();

    1037:             section->updateAddPoint(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
    1038:           }
    1039:         }
    1040:       }
    1041:     };
    1042:     class InsertionBinaryFusion {
    1043:     public:
    1044:       // Insert the overlapSection values into section along recvOverlap
    1045:       template<typename OverlapSection, typename RecvOverlap, typename Section>
    1046:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
    1047:         typedef typename Section::point_type        point_type;
    1048:         typedef typename OverlapSection::point_type overlap_point_type;
    1049: #if 0
    1050:         const Obj<typename RecvOverlap::baseSequence>      rPoints = recvOverlap->base();
    1051:         const typename RecvOverlap::baseSequence::iterator rEnd    = rPoints->end();

    1053:         for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
    1054:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
    1055:           const point_type&                              localPoint  = *p_iter;
    1056:           const int                                      rank        = *points->begin();
    1057:           const point_type&                              remotePoint = points->begin().color();

    1059:           if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
    1060:             if (!section->hasPoint(localPoint)) {
    1061:               std::cout <<"["<<section->commRank()<<"]: Destination section does not have local point " << localPoint << " remote point " << remotePoint << " fiber dim " << overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)) << std::endl;
    1062:             }
    1063:             section->setFiberDimension(localPoint, overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)));
    1064:           }
    1065:         }
    1066:         if (rPoints->size()) {section->allocatePoint();}
    1067:         for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
    1068:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
    1069:           const point_type&                              localPoint  = *p_iter;
    1070:           const int                                      rank        = *points->begin();
    1071:           const point_type&                              remotePoint = points->begin().color();

    1073:           if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
    1074:             section->updatePoint(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
    1075:           }
    1076:         }
    1077: #else
    1078:         const typename RecvOverlap::capSequence::iterator rBegin = recvOverlap->capBegin();
    1079:         const typename RecvOverlap::capSequence::iterator rEnd   = recvOverlap->capEnd();
    1080:         bool                                              hasPoints = false;

    1082:         for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
    1083:           const int                                             rank    = *r_iter;
    1084:           const typename RecvOverlap::supportSequence::iterator pBegin  = recvOverlap->supportBegin(*r_iter);
    1085:           const typename RecvOverlap::supportSequence::iterator pEnd    = recvOverlap->supportEnd(*r_iter);

    1087:           for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
    1088:             const point_type&                              localPoint  = *p_iter;
    1089:             const point_type&                              remotePoint = p_iter.color();

    1091:             if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
    1092:               if (!section->hasPoint(localPoint)) {
    1093:                 std::cout <<"["<<section->commRank()<<"]: Destination section does not have local point " << localPoint << " remote point " << remotePoint << " fiber dim " << overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)) << std::endl;
    1094:               }
    1095:               section->setFiberDimension(localPoint, overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)));
    1096:             }
    1097:             hasPoints = true;
    1098:           }
    1099:         }
    1100:         if (hasPoints) {section->allocatePoint();}
    1101:         for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
    1102:           const int                                             rank    = *r_iter;
    1103:           const typename RecvOverlap::supportSequence::iterator pBegin  = recvOverlap->supportBegin(*r_iter);
    1104:           const typename RecvOverlap::supportSequence::iterator pEnd    = recvOverlap->supportEnd(*r_iter);

    1106:           for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
    1107:             const point_type& localPoint  = *p_iter;
    1108:             const point_type& remotePoint = p_iter.color();

    1110:             if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
    1111:               section->updatePoint(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
    1112:             }
    1113:           }
    1114:         }
    1115: #endif
    1116:       }
    1117:       // Specialize to ArrowSection
    1118:       template<typename OverlapSection, typename RecvOverlap>
    1119:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<UniformSection<MinimalArrow<int,int>,int> >& section) {
    1120:         typedef UniformSection<MinimalArrow<int,int>,int> Section;
    1121:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
    1122:         const typename RecvOverlap::traits::baseSequence::iterator rBegin  = rPoints->begin();
    1123:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

    1125:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rBegin; p_iter != rEnd; ++p_iter) {
    1126:           const Obj<typename RecvOverlap::coneSequence>& sources      = recvOverlap->cone(*p_iter);
    1127:           const typename RecvOverlap::target_type        localSource  = *p_iter;
    1128:           const typename RecvOverlap::target_type        remoteSource = sources->begin().color();

    1130:           for(typename RecvOverlap::traits::baseSequence::iterator q_iter = rBegin; q_iter != rEnd; ++q_iter) {
    1131:             const Obj<typename RecvOverlap::coneSequence>& targets      = recvOverlap->cone(*q_iter);
    1132:             const typename RecvOverlap::target_type        localTarget  = *q_iter;
    1133:             const typename RecvOverlap::target_type        remoteTarget = targets->begin().color();
    1134:             const typename Section::point_type             localPoint(localSource, localTarget);
    1135:             const typename OverlapSection::point_type      remotePoint(remoteSource, remoteTarget);

    1137:             if (overlapSection->hasPoint(remotePoint)) {section->setFiberDimension(localPoint, overlapSection->getFiberDimension(remotePoint));}
    1138:           }
    1139:         }
    1140:         if (rPoints->size()) {section->allocatePoint();}
    1141:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rBegin; p_iter != rEnd; ++p_iter) {
    1142:           const Obj<typename RecvOverlap::coneSequence>& sources      = recvOverlap->cone(*p_iter);
    1143:           const typename RecvOverlap::target_type        localSource  = *p_iter;
    1144:           const typename RecvOverlap::target_type        remoteSource = sources->begin().color();

    1146:           for(typename RecvOverlap::traits::baseSequence::iterator q_iter = rBegin; q_iter != rEnd; ++q_iter) {
    1147:             const Obj<typename RecvOverlap::coneSequence>& targets      = recvOverlap->cone(*q_iter);
    1148:             const typename RecvOverlap::target_type        localTarget  = *q_iter;
    1149:             const typename RecvOverlap::target_type        remoteTarget = targets->begin().color();
    1150:             const typename Section::point_type             localPoint(localSource, localTarget);
    1151:             const typename OverlapSection::point_type      remotePoint(remoteSource, remoteTarget);

    1153:             if (overlapSection->hasPoint(remotePoint)) {section->updatePoint(localPoint, overlapSection->restrictPoint(remotePoint));}
    1154:           }
    1155:         }
    1156:       }
    1157:       // Specialize to the Sieve
    1158:       template<typename OverlapSection, typename RecvOverlap, typename Renumbering, typename Point>
    1159:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, Renumbering& renumbering, const Obj<Sieve<Point,Point,int> >& sieve) {
    1160:         typedef typename OverlapSection::point_type overlap_point_type;
    1161:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
    1162:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

    1164:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
    1165:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
    1166:           const Point&                                   localPoint  = *p_iter;
    1167:           const int                                      rank        = *points->begin();
    1168:           const Point&                                   remotePoint = points->begin().color();
    1169:           const int                                      size        = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
    1170:           const typename OverlapSection::value_type     *values      = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
    1171:           int                                            c           = 0;

    1173:           sieve->clearCone(localPoint);
    1174:           for(int i = 0; i < size; ++i, ++c) {sieve->addCone(renumbering[values[i]], localPoint, c);}
    1175:         }
    1176:       }
    1177:       // Specialize to the ISieve
    1178:       template<typename OverlapSection, typename RecvOverlap, typename Renumbering, typename Point>
    1179:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, Renumbering& renumbering, const Obj<IFSieve<Point> >& sieve) {
    1180:         typedef typename OverlapSection::point_type overlap_point_type;
    1181: #if 0
    1182:         const Obj<typename RecvOverlap::baseSequence>      rPoints = recvOverlap->base();
    1183:         const typename RecvOverlap::baseSequence::iterator rEnd    = rPoints->end();
    1184:         int                                                maxSize = 0;

    1186:         for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
    1187:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
    1188:           const Point&                                   localPoint  = *p_iter;
    1189:           const int                                      rank        = *points->begin();
    1190:           const Point&                                   remotePoint = points->begin().color();
    1191:           const int                                      size        = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
    1192:           const typename OverlapSection::value_type     *values      = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));

    1194:           sieve->setConeSize(localPoint, size);
    1195:           ///for(int i = 0; i < size; ++i) {sieve->addSupportSize(renumbering[values[i]], 1);}
    1196:           for(int i = 0; i < size; ++i) {sieve->addSupportSize(renumbering[values[i].first], 1);}
    1197:           maxSize = std::max(maxSize, size);
    1198:         }
    1199:         sieve->allocate();
    1200:         ///typename OverlapSection::value_type *localValues = new typename OverlapSection::value_type[maxSize];
    1201:         typename OverlapSection::value_type::first_type  *localValues      = new typename OverlapSection::value_type::first_type[maxSize];
    1202:         typename OverlapSection::value_type::second_type *localOrientation = new typename OverlapSection::value_type::second_type[maxSize];

    1204:         for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
    1205:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
    1206:           const Point&                                   localPoint  = *p_iter;
    1207:           const int                                      rank        = *points->begin();
    1208:           const Point&                                   remotePoint = points->begin().color();
    1209:           const int                                      size        = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
    1210:           const typename OverlapSection::value_type     *values      = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));

    1212:           ///for(int i = 0; i < size; ++i) {localValues[i] = renumbering[values[i]];}
    1213:           for(int i = 0; i < size; ++i) {
    1214:             localValues[i]      = renumbering[values[i].first];
    1215:             localOrientation[i] = values[i].second;
    1216:           }
    1217:           sieve->setCone(localValues, localPoint);
    1218:           sieve->setConeOrientation(localOrientation, localPoint);
    1219:         }
    1220:         delete [] localValues;
    1221:         delete [] localOrientation;
    1222: #else
    1223:         const typename RecvOverlap::capSequence::iterator rBegin  = recvOverlap->capBegin();
    1224:         const typename RecvOverlap::capSequence::iterator rEnd    = recvOverlap->capEnd();
    1225:         int                                               maxSize = 0;

    1227:         for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
    1228:           const int                                             rank    = *r_iter;
    1229:           const typename RecvOverlap::supportSequence::iterator pBegin  = recvOverlap->supportBegin(*r_iter);
    1230:           const typename RecvOverlap::supportSequence::iterator pEnd    = recvOverlap->supportEnd(*r_iter);

    1232:           for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
    1233:             const Point&                               localPoint  = *p_iter;
    1234:             const Point&                               remotePoint = p_iter.color();
    1235:             const int                                  size        = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
    1236:             const typename OverlapSection::value_type *values      = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));

    1238:             sieve->setConeSize(localPoint, size);
    1239:             for(int i = 0; i < size; ++i) {sieve->addSupportSize(renumbering[values[i].first], 1);}
    1240:             maxSize = std::max(maxSize, size);
    1241:           }
    1242:         }
    1243:         sieve->allocate();
    1244:         typename OverlapSection::value_type::first_type  *localValues      = new typename OverlapSection::value_type::first_type[maxSize];
    1245:         typename OverlapSection::value_type::second_type *localOrientation = new typename OverlapSection::value_type::second_type[maxSize];

    1247:         for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
    1248:           const int                                             rank    = *r_iter;
    1249:           const typename RecvOverlap::supportSequence::iterator pBegin  = recvOverlap->supportBegin(*r_iter);
    1250:           const typename RecvOverlap::supportSequence::iterator pEnd    = recvOverlap->supportEnd(*r_iter);

    1252:           for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
    1253:             const Point&                               localPoint  = *p_iter;
    1254:             const Point&                               remotePoint = p_iter.color();
    1255:             const int                                  size        = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
    1256:             const typename OverlapSection::value_type *values      = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));

    1258:             for(int i = 0; i < size; ++i) {
    1259:               localValues[i]      = renumbering[values[i].first];
    1260:               localOrientation[i] = values[i].second;
    1261:             }
    1262:             sieve->setCone(localValues, localPoint);
    1263:             sieve->setConeOrientation(localOrientation, localPoint);
    1264:           }
    1265:         }
    1266:         delete [] localValues;
    1267:         delete [] localOrientation;
    1268: #endif
    1269:       }
    1270:       template<typename OverlapSection, typename RecvOverlap, typename Point>
    1271:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<IFSieve<Point> >& sieve) {
    1272:         typedef typename OverlapSection::point_type overlap_point_type;
    1273: #if 0
    1274:         const Obj<typename RecvOverlap::baseSequence>      rPoints = recvOverlap->base();
    1275:         const typename RecvOverlap::baseSequence::iterator rEnd    = rPoints->end();
    1276:         int                                                maxSize = 0;

    1278:         for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
    1279:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
    1280:           const Point&                                   localPoint  = *p_iter;
    1281:           const int                                      rank        = *points->begin();
    1282:           const Point&                                   remotePoint = points->begin().color();
    1283:           const int                                      size        = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
    1284:           const typename OverlapSection::value_type     *values      = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));

    1286:           sieve->setConeSize(localPoint, size);
    1287:           for(int i = 0; i < size; ++i) {sieve->addSupportSize(values[i].first, 1);}
    1288:           maxSize = std::max(maxSize, size);
    1289:         }
    1290:         sieve->allocate();
    1291:         typename OverlapSection::value_type::first_type  *localValues      = new typename OverlapSection::value_type::first_type[maxSize];
    1292:         typename OverlapSection::value_type::second_type *localOrientation = new typename OverlapSection::value_type::second_type[maxSize];

    1294:         for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
    1295:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
    1296:           const Point&                                   localPoint  = *p_iter;
    1297:           const int                                      rank        = *points->begin();
    1298:           const Point&                                   remotePoint = points->begin().color();
    1299:           const int                                      size        = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
    1300:           const typename OverlapSection::value_type     *values      = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));

    1302:           for(int i = 0; i < size; ++i) {
    1303:             localValues[i]      = values[i].first;
    1304:             localOrientation[i] = values[i].second;
    1305:           }
    1306:           sieve->setCone(localValues, localPoint);
    1307:           sieve->setConeOrientation(localOrientation, localPoint);
    1308:         }
    1309:         delete [] localValues;
    1310:         delete [] localOrientation;
    1311: #else
    1312:         const Obj<typename RecvOverlap::capSequence>      rRanks  = recvOverlap->cap();
    1313:         const typename RecvOverlap::capSequence::iterator rEnd    = rRanks->end();
    1314:         int                                               maxSize = 0;

    1316:         for(typename RecvOverlap::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
    1317:           const int                                             rank    = *r_iter;
    1318:           const Obj<typename RecvOverlap::supportSequence>&     points  = recvOverlap->support(*r_iter);
    1319:           const typename RecvOverlap::supportSequence::iterator pBegin  = points->begin();
    1320:           const typename RecvOverlap::supportSequence::iterator pEnd    = points->end();

    1322:           for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
    1323:             const Point&                               localPoint  = *p_iter;
    1324:             const Point&                               remotePoint = p_iter.color();
    1325:             const int                                  size        = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
    1326:             const typename OverlapSection::value_type *values      = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));

    1328:             sieve->setConeSize(localPoint, size);
    1329:             for(int i = 0; i < size; ++i) {sieve->addSupportSize(values[i].first, 1);}
    1330:             maxSize = std::max(maxSize, size);
    1331:           }
    1332:         }
    1333:         sieve->allocate();
    1334:         typename OverlapSection::value_type::first_type  *localValues      = new typename OverlapSection::value_type::first_type[maxSize];
    1335:         typename OverlapSection::value_type::second_type *localOrientation = new typename OverlapSection::value_type::second_type[maxSize];

    1337:         for(typename RecvOverlap::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
    1338:           const int                                             rank    = *r_iter;
    1339:           const Obj<typename RecvOverlap::supportSequence>&     points  = recvOverlap->support(*r_iter);
    1340:           const typename RecvOverlap::supportSequence::iterator pBegin  = points->begin();
    1341:           const typename RecvOverlap::supportSequence::iterator pEnd    = points->end();

    1343:           for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
    1344:             const Point&                               localPoint  = *p_iter;
    1345:             const Point&                               remotePoint = p_iter.color();
    1346:             const int                                  size        = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
    1347:             const typename OverlapSection::value_type *values      = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));

    1349:             for(int i = 0; i < size; ++i) {
    1350:               localValues[i]      = values[i].first;
    1351:               localOrientation[i] = values[i].second;
    1352:             }
    1353:             sieve->setCone(localValues, localPoint);
    1354:             sieve->setConeOrientation(localOrientation, localPoint);
    1355:           }
    1356:         }
    1357:         delete [] localValues;
    1358:         delete [] localOrientation;
    1359: #endif
    1360:       }
    1361:       // Generic
    1362:       template<typename OverlapSection, typename RecvOverlap, typename Section, typename Bundle>
    1363:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section, const Obj<Bundle>& bundle) {
    1364:         typedef typename OverlapSection::point_type overlap_point_type;
    1365:         const Obj<typename RecvOverlap::baseSequence>      rPoints = recvOverlap->base();
    1366:         const typename RecvOverlap::baseSequence::iterator rEnd    = rPoints->end();

    1368:         for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
    1369:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
    1370:           const typename Section::point_type&            localPoint  = *p_iter;
    1371:           const int                                      rank        = *points->begin();
    1372:           const typename OverlapSection::point_type&     remotePoint = points->begin().color();

    1374:           section->setFiberDimension(localPoint, overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)));
    1375:         }
    1376:         bundle->allocate(section);
    1377:         for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
    1378:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
    1379:           const typename Section::point_type&            localPoint  = *p_iter;
    1380:           const int                                      rank        = *points->begin();
    1381:           const typename OverlapSection::point_type&     remotePoint = points->begin().color();

    1383:           section->update(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
    1384:         }
    1385:       }
    1386:     };
    1387:     class InterpolateMultipleFusion {
    1388:     public:
    1389:       // Interpolate the overlapSection values into section along recvOverlap
    1390:       template<typename OverlapSection, typename RecvOverlap, typename Section>
    1391:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
    1392:         typedef typename Section::point_type        point_type;
    1393:         typedef typename Section::value_type        value_type;
    1394:         typedef typename OverlapSection::point_type overlap_point_type;
    1395:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints     = recvOverlap->base();
    1396:         const typename RecvOverlap::traits::baseSequence::iterator rEnd        = rPoints->end();
    1397:         int                                                        maxFiberDim = -1;

    1399:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
    1400:           const Obj<typename RecvOverlap::coneSequence>&     points     = recvOverlap->cone(*p_iter);
    1401:           const typename RecvOverlap::coneSequence::iterator rpEnd      = points->end();
    1402:           const point_type&                                  localPoint = *p_iter;
    1403:           bool                                               inOverlap  = false;
    1404:           int                                                fiberDim   = -1;

    1406:           for(typename RecvOverlap::coneSequence::iterator rp_iter = points->begin(); rp_iter != rpEnd; ++rp_iter) {
    1407:             const int         rank        = *rp_iter;
    1408:             const point_type& remotePoint = rp_iter.color();

    1410:             if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
    1411:               inOverlap = true;
    1412:               fiberDim  = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
    1413:               break;
    1414:             }
    1415:           }
    1416:           if (inOverlap) {
    1417:             if (!section->hasPoint(localPoint)) {
    1418:               std::cout <<"["<<section->commRank()<<"]: Destination section does not have local point " << localPoint << " remote point " << (points->begin().color()) << " fiber dim " << fiberDim << std::endl;
    1419:             }
    1420:             section->setFiberDimension(localPoint, fiberDim);
    1421:             maxFiberDim = std::max(fiberDim, maxFiberDim);
    1422:           }
    1423:         }
    1424:         if (rPoints->size()) {section->allocatePoint();}
    1425:         value_type *interpolant = new value_type[maxFiberDim];

    1427:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
    1428:           const Obj<typename RecvOverlap::coneSequence>&     points     = recvOverlap->cone(*p_iter);
    1429:           const typename RecvOverlap::coneSequence::iterator rpEnd      = points->end();
    1430:           const point_type&                                  localPoint = *p_iter;
    1431:           bool                                               inOverlap  = false;
    1432:           int                                                numArgs    = 0;

    1434:           for(int d = 0; d < maxFiberDim; ++d) {interpolant[d] = 0.0;}
    1435:           for(typename RecvOverlap::coneSequence::iterator rp_iter = points->begin(); rp_iter != rpEnd; ++rp_iter) {
    1436:             const int         rank        = *rp_iter;
    1437:             const point_type& remotePoint = rp_iter.color();
    1438:             const overlap_point_type opoint(rank, remotePoint);

    1440:             if (overlapSection->hasPoint(opoint)) {
    1441:               const int         fiberDim = overlapSection->getFiberDimension(opoint);
    1442:               const value_type *values   = overlapSection->restrictPoint(opoint);

    1444:               // TODO: Include interpolation weights (stored in overlap)
    1445:               for(int d = 0; d < fiberDim; ++d) {
    1446:                 interpolant[d] += values[d];
    1447:               }
    1448:               inOverlap = true;
    1449:               ++numArgs;
    1450:             }
    1451:           }
    1452:           if (inOverlap) {
    1453:             for(int d = 0; d < maxFiberDim; ++d) {interpolant[d] /= numArgs;}
    1454:             section->updatePoint(localPoint, interpolant);
    1455:           }
    1456:         }
    1457:         delete [] interpolant;
    1458:       }
    1459:     };
    1460:   }
    1461: }

    1463: #endif