Actual source code: ParallelMapping.hh

  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:     }
344:   };
345:   namespace Pullback {
346:     class SimpleCopy {
347:     public:
348:       // Copy the overlap section to the related processes
349:       //   This version is for Constant sections, meaning the same, single value over all points
350:       template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
351:       static void copyConstant(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
352:         MPIMover<char>                             pMover(sendSection->comm(), sendSection->debug());
353:         MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), sendSection->debug());
354:         std::map<int, char *>                      sendPoints;
355:         std::map<int, char *>                      recvPoints;
356:         typename SendSection::alloc_type::template rebind<char>::other sendAllocator;
357:         typename RecvSection::alloc_type::template rebind<char>::other recvAllocator;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

525:         typename SendSection::point_type min = -1;
526:         typename SendSection::point_type max = -1;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1057:           if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1058:             if (!section->hasPoint(localPoint)) {
1059:               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;
1060:             }
1061:             section->setFiberDimension(localPoint, overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)));
1062:           }
1063:         }
1064:         if (rPoints->size()) {section->allocatePoint();}
1065:         for(typename RecvOverlap::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1066:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
1067:           const point_type&                              localPoint  = *p_iter;
1068:           const int                                      rank        = *points->begin();
1069:           const point_type&                              remotePoint = points->begin().color();

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

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

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

1089:             if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1090:               if (!section->hasPoint(localPoint)) {
1091:                 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;
1092:               }
1093:               section->setFiberDimension(localPoint, overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)));
1094:             }
1095:             hasPoints = true;
1096:           }
1097:         }
1098:         if (hasPoints) {section->allocatePoint();}
1099:         for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
1100:           const int                                             rank    = *r_iter;
1101:           const typename RecvOverlap::supportSequence::iterator pBegin  = recvOverlap->supportBegin(*r_iter);
1102:           const typename RecvOverlap::supportSequence::iterator pEnd    = recvOverlap->supportEnd(*r_iter);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1461: #endif