Actual source code: ParallelMapping.hh

  1: #ifndef included_ALE_ParallelMapping_hh
  2: #define included_ALE_ParallelMapping_hh

  4: #ifndef  included_ALE_IField_hh
  5: #include <IField.hh>
  6: #endif

  8: #ifndef  included_ALE_Sections_hh
  9: #include <Sections.hh>
 10: #endif

 12: #include <functional>


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

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

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

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

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

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

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

156:   // TODO: Check MPI return values and status of Waits
157:   template<typename Value_>
158:   class MPIMover : public ALE::ParallelObject {
159:   public:
160:     typedef Value_                                 value_type;
161:     typedef size_t                                 num_type;
162:     typedef std::pair<num_type, const value_type*> move_type;
163:     typedef std::map<int, move_type>               moves_type;
164:     typedef std::vector<MPI_Request>               requests_type;
165:   protected:
166:     bool          _createdType;
167:     int           _tag;
168:     MPI_Datatype  _datatype;
169:     moves_type    _sends;
170:     moves_type    _recvs;
171:     requests_type _requests;
172:   public:
173:     MPIMover(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug), _createdType(0) {
174:       this->_tag      = this->getNewTag();
175:       this->_datatype = this->getMPIDatatype();
176:     };
177:     MPIMover(MPI_Comm comm, const int tag, const int debug) : ParallelObject(comm, debug), _createdType(0), _tag(tag) {
178:       this->_datatype = this->getMPIDatatype();
179:     };
180:     MPIMover(MPI_Comm comm, const MPI_Datatype datatype, const int tag, const int debug) : ParallelObject(comm, debug), _createdType(0) {
181:       if (tag == MPI_UNDEFINED) {
182:         this->_tag      = this->getNewTag();
183:       } else {
184:         this->_tag      = tag;
185:       }
186:       if (datatype == MPI_DATATYPE_NULL) {
187:         this->_datatype = this->getMPIDatatype();
188:       } else {
189:         this->_datatype = datatype;
190:       }
191:     };
192:     ~MPIMover() {
193:       if (_createdType) {
194:         MPI_Type_free(&this->_datatype);
195:       }
196:     };
197:   protected:
198:     // TODO: Can rewrite this with template specialization?
199:     MPI_Datatype getMPIDatatype() {
200:       if (sizeof(value_type) == 1) {
201:         return MPI_BYTE;
202:       } else if (sizeof(value_type) == 2) {
203:         return MPI_SHORT;
204:       } else if (sizeof(value_type) == 4) {
205:         return MPI_INT;
206:       } else if (sizeof(value_type) == 8) {
207:         return MPI_DOUBLE;
208:       } else if (sizeof(value_type) == 28) {
209:         int          blen[2];
210:         MPI_Aint     indices[2];
211:         MPI_Datatype oldtypes[2], newtype;
212:         blen[0] = 1; indices[0] = 0;           oldtypes[0] = MPI_INT;
213:         blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
214:         MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
215:         MPI_Type_commit(&newtype);
216:         this->_createdType = true;
217:         return newtype;
218:       } else if (sizeof(value_type) == 32) {
219:         int          blen[2];
220:         MPI_Aint     indices[2];
221:         MPI_Datatype oldtypes[2], newtype;
222:         blen[0] = 1; indices[0] = 0;           oldtypes[0] = MPI_DOUBLE;
223:         blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
224:         MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
225:         MPI_Type_commit(&newtype);
226:         this->_createdType = true;
227:         return newtype;
228:       }
229:       throw ALE::Exception("Cannot determine MPI type for value type");
230:     };
231:     int getNewTag() const {
232:       static int tagKeyval = MPI_KEYVAL_INVALID;
233:       int *tagvalp = NULL, *maxval, flg;

235:       if (tagKeyval == MPI_KEYVAL_INVALID) {
236:         tagvalp = (int *) malloc(sizeof(int));
237:         MPI_Keyval_create(MPI_NULL_COPY_FN, Mesh_DelTag, &tagKeyval, (void *) NULL);
238:         MPI_Attr_put(this->comm(), tagKeyval, tagvalp);
239:         tagvalp[0] = 0;
240:       }
241:       MPI_Attr_get(this->comm(), tagKeyval, (void **) &tagvalp, &flg);
242:       if (tagvalp[0] < 1) {
243:         MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, (void **) &maxval, &flg);
244:         tagvalp[0] = *maxval - 128; // hope that any still active tags were issued right at the beginning of the run
245:       }
246:       if (this->debug()) {
247:         std::cout << "[" << this->commRank() << "]Got new tag " << tagvalp[0] << std::endl;
248:       }
249:       return tagvalp[0]--;
250:     };
251:     void constructRequests() {
252:       this->_requests.clear();

254:       for(typename moves_type::const_iterator s_iter = this->_sends.begin(); s_iter != this->_sends.end(); ++s_iter) {
255:         const int   rank = s_iter->first;
256:         const int   num  = s_iter->second.first;
257:         void       *data = (void *) s_iter->second.second;
258:         MPI_Request request;

260:         if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Sending data (" << num << ") to " << rank << " tag " << this->_tag << std::endl;}
261:         MPI_Send_init(data, num, this->_datatype, rank, this->_tag, this->comm(), &request);
262:         this->_requests.push_back(request);
263: #if defined(PETSC_USE_LOG)
264:         // PETSc logging
265:         isend_ct++;
266:         TypeSize(&isend_len, num, this->_datatype);
267: #endif
268:       }
269:       for(typename moves_type::const_iterator r_iter = this->_recvs.begin(); r_iter != this->_recvs.end(); ++r_iter) {
270:         const int   rank = r_iter->first;
271:         const int   num  = r_iter->second.first;
272:         void       *data = (void *) r_iter->second.second;
273:         MPI_Request request;

275:         if (this->_debug) {std::cout <<"["<<this->commRank()<<"] Receiving data (" << num << ") from " << rank << " tag " << this->_tag << std::endl;}
276:         MPI_Recv_init(data, num, this->_datatype, rank, this->_tag, this->comm(), &request);
277:         this->_requests.push_back(request);
278: #if defined(PETSC_USE_LOG)
279:         // PETSc logging
280:         irecv_ct++;
281:         TypeSize(&irecv_len, num, this->_datatype);
282: #endif
283:       }
284:     };
285:   public:
286:     void send(const int rank, const int num, const value_type *data) {
287:       this->_sends[rank] = move_type(num, data);
288:     };
289:     void recv(const int rank, const int num, const value_type *data) {
290:       this->_recvs[rank] = move_type(num, data);
291:     };
292:     void start() {
293:       this->constructRequests();
294:       for(typename requests_type::const_iterator r_iter = this->_requests.begin(); r_iter != this->_requests.end(); ++r_iter) {
295:         MPI_Request request = *r_iter;

297:         MPI_Start(&request);
298:       }
299:     };
300:     void end() {
301:       MPI_Status status;

303:       for(typename requests_type::const_iterator r_iter = this->_requests.begin(); r_iter != this->_requests.end(); ++r_iter) {
304:         MPI_Request request = *r_iter;

306:         MPI_Wait(&request, &status);
307:       }
308:     };
309:   };
310:   template<typename Alloc_ = malloc_allocator<int> >
311:   class OverlapBuilder {
312:   public:
313:     typedef Alloc_ alloc_type;
314:   protected:
315:     template<typename T>
316:     struct lessPair: public std::binary_function<std::pair<T,T>, std::pair<T,T>, bool> {
317:       bool operator()(const std::pair<T,T>& x, const std::pair<T,T>& y) const {
318:         return x.first < y.first;
319:       }
320:     };
321:     template<typename T>
322:     struct mergePair: public std::binary_function<std::pair<T,T>, std::pair<T,T>, bool> {
323:       std::pair<T,std::pair<T,T> > operator()(const std::pair<T,T>& x, const std::pair<T,T>& y) const {
324:         return std::pair<T,std::pair<T,T> >(x.first, std::pair<T,T>(x.second, y.second));
325:       }
326:     };
327:     template<typename _InputIterator1, typename _InputIterator2, typename _OutputIterator, typename _Compare, typename _Merge>
328:     static _OutputIterator set_intersection_merge(_InputIterator1 __first1, _InputIterator1 __last1,
329:                                            _InputIterator2 __first2, _InputIterator2 __last2,
330:                                            _OutputIterator __result, _Compare __comp, _Merge __merge)
331:     {
332:       while(__first1 != __last1 && __first2 != __last2) {
333:         if (__comp(*__first1, *__first2))
334:           ++__first1;
335:         else if (__comp(*__first2, *__first1))
336:           ++__first2;
337:         else
338:         {
339:           *__result = __merge(*__first1, *__first2);
340:           ++__first1;
341:           ++__first2;
342:           ++__result;
343:         }
344:       }
345:       return __result;
346:     };
347:   public:
348:     template<typename Sequence, typename Renumbering, typename SendOverlap, typename RecvOverlap>
349:     static void constructOverlap(const Sequence& points, Renumbering& renumbering, const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap) {
350:       typedef typename SendOverlap::source_type point_type;
351:       typedef std::pair<point_type,point_type>  pointPair;
352:       typedef std::pair<point_type,pointPair>   pointTriple;
353:       alloc_type allocator;
354:       typename alloc_type::template rebind<point_type>::other point_allocator;
355:       typename alloc_type::template rebind<pointPair>::other  pointPair_allocator;
356:       const MPI_Comm comm     = sendOverlap->comm();
357:       const int      commSize = sendOverlap->commSize();
358:       const int      commRank = sendOverlap->commRank();
359:       point_type    *sendBuf  = point_allocator.allocate(points.size()*2);
360:       for(size_t i = 0; i < points.size()*2; ++i) {point_allocator.construct(sendBuf+i, point_type());}
361:       int            size     = 0;
362:       const int      debug    = sendOverlap->debug();
363:       for(typename Sequence::const_iterator l_iter = points.begin(); l_iter != points.end(); ++l_iter) {
364:         if (debug) {std::cout << "["<<commRank<<"]Send point["<<size<<"]: " << *l_iter << " " << renumbering[*l_iter] << std::endl;}
365:         sendBuf[size++] = *l_iter;
366:         sendBuf[size++] = renumbering[*l_iter];
367:       }
368:       int *sizes = allocator.allocate(commSize*3+2); // [size]   The number of points coming from each process
369:       for(int i = 0; i < commSize*3+2; ++i) {allocator.construct(sizes+i, 0);}
370:       int *offsets = sizes+commSize;                 // [size+1] Prefix sums for sizes
371:       int *oldOffs = offsets+commSize+1;             // [size+1] Temporary storage
372:       pointPair  *remotePoints = NULL;               // The points from each process
373:       int        *remoteRanks  = NULL;               // The rank and number of overlap points of each process that overlaps another
374:       int         numRemotePoints = 0;
375:       int         numRemoteRanks  = 0;

377:       // Change to Allgather() for the correct binning algorithm
378:       MPI_Gather(&size, 1, MPI_INT, sizes, 1, MPI_INT, 0, comm);
379:       if (commRank == 0) {
380:         offsets[0] = 0;
381:         for(int p = 1; p <= commSize; p++) {
382:           offsets[p] = offsets[p-1] + sizes[p-1];
383:         }
384:         numRemotePoints = offsets[commSize];
385:         remotePoints    = pointPair_allocator.allocate(numRemotePoints/2);
386:         for(int i = 0; i < numRemotePoints/2; ++i) {pointPair_allocator.construct(remotePoints+i, pointPair());}
387:       }
388:       MPI_Gatherv(sendBuf, size, MPI_INT, remotePoints, sizes, offsets, MPI_INT, 0, comm);
389:       for(size_t i = 0; i < points.size(); ++i) {point_allocator.destroy(sendBuf+i);}
390:       point_allocator.deallocate(sendBuf, points.size());
391:       std::map<int, std::map<int, std::set<pointTriple> > > overlapInfo; // Maps (p,q) to their set of overlap points

393:       if (commRank == 0) {
394:         for(int p = 0; p <= commSize; p++) {
395:           offsets[p] /= 2;
396:         }
397:         for(int p = 0; p < commSize; p++) {
398:           std::sort(&remotePoints[offsets[p]], &remotePoints[offsets[p+1]], lessPair<point_type>());
399:         }
400:         for(int p = 0; p <= commSize; p++) {
401:           oldOffs[p] = offsets[p];
402:         }
403:         for(int p = 0; p < commSize; p++) {
404:           for(int q = 0; q < commSize; q++) {
405:             if (p == q) continue;
406:             set_intersection_merge(&remotePoints[oldOffs[p]], &remotePoints[oldOffs[p+1]],
407:                                    &remotePoints[oldOffs[q]], &remotePoints[oldOffs[q+1]],
408:                                    std::insert_iterator<std::set<pointTriple> >(overlapInfo[p][q], overlapInfo[p][q].begin()),
409:                                    lessPair<point_type>(), mergePair<point_type>());
410:           }
411:           sizes[p]     = overlapInfo[p].size()*2;
412:           offsets[p+1] = offsets[p] + sizes[p];
413:         }
414:         numRemoteRanks = offsets[commSize];
415:         remoteRanks    = allocator.allocate(numRemoteRanks);
416:         for(int i = 0; i < numRemoteRanks; ++i) {allocator.construct(remoteRanks+i, 0);}
417:         int     k = 0;
418:         for(int p = 0; p < commSize; p++) {
419:           for(typename std::map<int, std::set<pointTriple> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
420:             remoteRanks[k*2]   = r_iter->first;
421:             remoteRanks[k*2+1] = r_iter->second.size();
422:             k++;
423:           }
424:         }
425:       }
426:       int numOverlaps;                          // The number of processes overlapping this process
427:       MPI_Scatter(sizes, 1, MPI_INT, &numOverlaps, 1, MPI_INT, 0, comm);
428:       int *overlapRanks = allocator.allocate(numOverlaps); // The rank and overlap size for each overlapping process
429:       for(int i = 0; i < numOverlaps; ++i) {allocator.construct(overlapRanks+i, 0);}
430:       MPI_Scatterv(remoteRanks, sizes, offsets, MPI_INT, overlapRanks, numOverlaps, MPI_INT, 0, comm);
431:       point_type *sendPoints    = NULL;         // The points to send to each process
432:       int         numSendPoints = 0;
433:       if (commRank == 0) {
434:         for(int p = 0, k = 0; p < commSize; p++) {
435:           sizes[p] = 0;
436:           for(int r = 0; r < (int) overlapInfo[p].size(); r++) {
437:             sizes[p] += remoteRanks[k*2+1]*2;
438:             k++;
439:           }
440:           offsets[p+1] = offsets[p] + sizes[p];
441:         }
442:         numSendPoints = offsets[commSize];
443:         sendPoints    = point_allocator.allocate(numSendPoints*2);
444:         for(int i = 0; i < numSendPoints*2; ++i) {point_allocator.construct(sendPoints+i, point_type());}
445:         for(int p = 0, k = 0; p < commSize; p++) {
446:           for(typename std::map<int, std::set<pointTriple> >::const_iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
447:             int rank = r_iter->first;
448:             for(typename std::set<pointTriple>::const_iterator p_iter = (overlapInfo[p][rank]).begin(); p_iter != (overlapInfo[p][rank]).end(); ++p_iter) {
449:               sendPoints[k++] = p_iter->first;
450:               sendPoints[k++] = p_iter->second.second;
451:               if (debug) {std::cout << "["<<commRank<<"]Sending points " << p_iter->first << " " << p_iter->second.second << " to rank " << rank << std::endl;}
452:             }
453:           }
454:         }
455:       }
456:       int numOverlapPoints = 0;
457:       for(int r = 0; r < numOverlaps/2; r++) {
458:         numOverlapPoints += overlapRanks[r*2+1];
459:       }
460:       point_type *overlapPoints = point_allocator.allocate(numOverlapPoints*2);
461:       for(int i = 0; i < numOverlapPoints*2; ++i) {point_allocator.construct(overlapPoints+i, point_type());}
462:       MPI_Scatterv(sendPoints, sizes, offsets, MPI_INT, overlapPoints, numOverlapPoints*2, MPI_INT, 0, comm);

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

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

471:           if (debug) {std::cout << "["<<commRank<<"]Matched up remote point " << remotePoint << "("<<point<<") to local " << renumbering[point] << std::endl;}
472:           sendOverlap->addArrow(renumbering[point], rank, remotePoint);
473:           recvOverlap->addArrow(rank, renumbering[point], remotePoint);
474:         }
475:       }

477:       for(int i = 0; i < numOverlapPoints; ++i) {point_allocator.destroy(overlapPoints+i);}
478:       point_allocator.deallocate(overlapPoints, numOverlapPoints);
479:       for(int i = 0; i < numOverlaps; ++i) {allocator.destroy(overlapRanks+i);}
480:       allocator.deallocate(overlapRanks, numOverlaps);
481:       for(int i = 0; i < commSize*3+2; ++i) {allocator.destroy(sizes+i);}
482:       allocator.deallocate(sizes, commSize*3+2);
483:       if (commRank == 0) {
484:         for(int i = 0; i < numRemoteRanks; ++i) {allocator.destroy(remoteRanks+i);}
485:         allocator.deallocate(remoteRanks, numRemoteRanks);
486:         for(int i = 0; i < numRemotePoints; ++i) {pointPair_allocator.destroy(remotePoints+i);}
487:         pointPair_allocator.deallocate(remotePoints, numRemotePoints);
488:         for(int i = 0; i < numSendPoints; ++i) {point_allocator.destroy(sendPoints+i);}
489:         point_allocator.deallocate(sendPoints, numSendPoints);
490:       }
491:     };
492:   };
493:   namespace Pullback {
494:     class SimpleCopy {
495:     public:
496:       // Copy the overlap section to the related processes
497:       //   This version is for Constant sections, meaning the same, single value over all points
498:       template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
499:       static void copyConstant(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
500:         MPIMover<char>                             pMover(sendSection->comm(), sendSection->debug());
501:         MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), sendSection->debug());
502:         std::map<int, char *>                      sendPoints;
503:         std::map<int, char *>                      recvPoints;
504:         typename SendSection::alloc_type::template rebind<char>::other sendAllocator;
505:         typename RecvSection::alloc_type::template rebind<char>::other recvAllocator;

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

511:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
512:           const Obj<typename SendOverlap::coneSequence>&     points = sendOverlap->cone(*r_iter);
513:           const int                                          pSize  = points->size();
514:           const typename SendOverlap::coneSequence::iterator pEnd   = points->end();
515:           char                                              *v      = sendAllocator.allocate(points->size());
516:           int                                                k      = 0;

518:           for(int i = 0; i < pSize; ++i) {sendAllocator.construct(v+i, 0);}
519:           for(typename SendOverlap::coneSequence::iterator p_iter = points->begin(); p_iter != pEnd; ++p_iter, ++k) {
520:             v[k] = (char) sendSection->hasPoint(*p_iter);
521:           }
522:           sendPoints[*r_iter] = v;
523:           pMover.send(*r_iter, pSize, sendPoints[*r_iter]);
524:           vMover.send(*r_iter, 2, sValues);
525:         }
526:         const Obj<typename RecvOverlap::traits::capSequence>      rRanks  = recvOverlap->cap();
527:         const typename RecvOverlap::traits::capSequence::iterator rEnd    = rRanks->end();
528:         const typename RecvSection::value_type                   *rValues = recvSection->restrictSpace();

530:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
531:           const Obj<typename RecvOverlap::traits::supportSequence>& points = recvOverlap->support(*r_iter);
532:           const int                                                 pSize  = points->size();
533:           char                                                     *v      = recvAllocator.allocate(pSize);

535:           for(int i = 0; i < pSize; ++i) {recvAllocator.construct(v+i, 0);}
536:           recvPoints[*r_iter] = v;
537:           pMover.recv(*r_iter, pSize, recvPoints[*r_iter]);
538:           vMover.recv(*r_iter, 2, rValues);
539:         }
540:         pMover.start();
541:         pMover.end();
542:         vMover.start();
543:         vMover.end();
544:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
545:           const Obj<typename RecvOverlap::traits::supportSequence>&     points = recvOverlap->support(*r_iter);
546:           const typename RecvOverlap::traits::supportSequence::iterator pEnd   = points->end();
547:           const char                                                   *v      = recvPoints[*r_iter];
548:           int                                                           k      = 0;

550:           for(typename RecvOverlap::traits::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter, ++k) {
551:             if (v[k]) {recvSection->addPoint(typename RecvSection::point_type(*r_iter, s_iter.color()));}
552:           }
553:         }
554:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
555:           sendAllocator.deallocate(sendPoints[*r_iter], sendOverlap->cone(*r_iter)->size());
556:         }
557:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
558:           recvAllocator.deallocate(recvPoints[*r_iter], recvOverlap->support(*r_iter)->size());
559:         }
560:       };
561:       // Specialize to ArrowSection
562:       template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
563:       static void copyConstantArrow(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
564:         MPIMover<char>                             pMover(sendSection->comm(), sendSection->debug());
565:         MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), sendSection->debug());
566:         std::map<int, char *>                      sendPoints;
567:         std::map<int, char *>                      recvPoints;
568:         typename SendSection::alloc_type::template rebind<char>::other sendAllocator;
569:         typename RecvSection::alloc_type::template rebind<char>::other recvAllocator;

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

575:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
576:           const Obj<typename SendOverlap::coneSequence>&     points = sendOverlap->cone(*r_iter);
577:           const int                                          pSize  = points->size();
578:           const typename SendOverlap::coneSequence::iterator pBegin = points->begin();
579:           const typename SendOverlap::coneSequence::iterator pEnd   = points->end();
580:           char                                              *v      = sendAllocator.allocate(pSize*pSize);
581:           int                                                k      = 0;

583:           for(size_t i = 0; i < pSize*pSize; ++i) {sendAllocator.construct(v+i, 0);}
584:           for(typename SendOverlap::coneSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
585:             for(typename SendOverlap::coneSequence::iterator q_iter = pBegin; q_iter != pEnd; ++q_iter, ++k) {
586:               v[k] = (char) sendSection->hasPoint(typename SendSection::point_type(*p_iter,*q_iter));
587:             }
588:           }
589:           sendPoints[*r_iter] = v;
590:           pMover.send(*r_iter, pSize*pSize, sendPoints[*r_iter]);
591:           vMover.send(*r_iter, 2, sValues);
592:         }
593:         const Obj<typename RecvOverlap::traits::capSequence>      rRanks  = recvOverlap->cap();
594:         const typename RecvOverlap::traits::capSequence::iterator rEnd    = rRanks->end();
595:         const typename RecvSection::value_type                   *rValues = recvSection->restrictSpace();

597:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
598:           const Obj<typename RecvOverlap::traits::supportSequence>& points = recvOverlap->support(*r_iter);
599:           const int                                                 pSize  = points->size();
600:           char                                                     *v      = recvAllocator.allocate(points->size()*points->size());

602:           for(size_t i = 0; i < pSize*pSize; ++i) {recvAllocator.construct(v+i, 0);}
603:           recvPoints[*r_iter] = v;
604:           pMover.recv(*r_iter, pSize*pSize, recvPoints[*r_iter]);
605:           vMover.recv(*r_iter, 2, rValues);
606:         }
607:         pMover.start();
608:         pMover.end();
609:         vMover.start();
610:         vMover.end();
611:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
612:           const Obj<typename RecvOverlap::traits::supportSequence>&     points = recvOverlap->support(*r_iter);
613:           const typename RecvOverlap::traits::supportSequence::iterator pBegin = points->begin();
614:           const typename RecvOverlap::traits::supportSequence::iterator pEnd   = points->end();
615:           const char                                                   *v      = recvPoints[*r_iter];
616:           int                                                           k      = 0;

618:           for(typename RecvOverlap::traits::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
619:             for(typename RecvOverlap::traits::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter, ++k) {
620:               if (v[k]) {recvSection->addPoint(typename RecvSection::point_type(s_iter.color(),t_iter.color()));}
621:             }
622:           }
623:         }
624:       };
625:       // Copy the overlap section to the related processes
626:       //   This version is for IConstant sections, meaning the same, single value over all points
627:       template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
628:       static void copyIConstant(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
629:         MPIMover<typename SendSection::point_type> pMover(sendSection->comm(), sendSection->debug());
630:         MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), sendSection->debug());
631:         std::map<int, typename SendSection::point_type *> sendPoints;
632:         std::map<int, typename SendSection::point_type *> recvPoints;
633:         typename SendSection::alloc_type::template rebind<typename SendSection::point_type>::other sendAllocator;
634:         typename RecvSection::alloc_type::template rebind<typename SendSection::point_type>::other recvAllocator;

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

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

643:           for(size_t i = 0; i < 2; ++i) {sendAllocator.construct(v+i, 0);}
644:           v[0] = sendSection->getChart().min();
645:           v[1] = sendSection->getChart().max();
646:           sendPoints[*r_iter] = v;
647:           pMover.send(*r_iter, 2, sendPoints[*r_iter]);
648:           vMover.send(*r_iter, 2, sValues);
649:           ///std::cout << "["<<sendOverlap->commRank()<<"]Sending chart (" << v[0] << ", " << v[1] << ") with values (" << sValues[0] << ", " << sValues[1] << ") to process " << *r_iter << std::endl;
650:         }
651:         const Obj<typename RecvOverlap::traits::capSequence>      rRanks  = recvOverlap->cap();
652:         const typename RecvOverlap::traits::capSequence::iterator rEnd    = rRanks->end();
653:         const typename RecvSection::value_type                   *rValues = recvSection->restrictSpace();

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

658:           for(size_t i = 0; i < 2; ++i) {recvAllocator.construct(v+i, 0);}
659:           recvPoints[*r_iter] = v;
660:           pMover.recv(*r_iter, 2, recvPoints[*r_iter]);
661:           vMover.recv(*r_iter, 2, rValues);
662:         }
663:         pMover.start();
664:         pMover.end();
665:         vMover.start();
666:         vMover.end();

668:         typename SendSection::point_type min = -1;
669:         typename SendSection::point_type max = -1;

671:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
672:           const typename RecvSection::point_type *v = recvPoints[*r_iter];
673:           typename SendSection::point_type        newMin = v[0];
674:           typename SendSection::point_type        newMax = v[1]-1;
675:           //int                                     pSize  = 0;

677:           ///std::cout << "["<<recvOverlap->commRank()<<"]Received chart (" << v[0] << ", " << v[1] << ") from process " << *r_iter << std::endl;
678: #if 0
679:           // Translate to local numbering
680:           if (recvOverlap->support(*r_iter)->size()) {
681:             while(!pSize) {
682:               const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter, newMin);
683:               pSize = points->size();
684:               if (pSize) {
685:                 newMin = *points->begin();
686:               } else {
687:                 newMin++;
688:               }
689:             }
690:             pSize  = 0;
691:             while(!pSize) {
692:               const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter, newMax);
693:               pSize = points->size();
694:               if (pSize) {
695:                 newMax = *points->begin();
696:               } else {
697:                 newMax--;
698:               }
699:             }
700:           }
701:           std::cout << "["<<recvOverlap->commRank()<<"]Translated to chart (" << newMin << ", " << newMax+1 << ") from process " << *r_iter << std::endl;
702: #endif
703:           // Update chart
704:           if (min < 0) {
705:             min = newMin;
706:             max = newMax+1;
707:           } else {
708:             min = std::min(min, newMin);
709:             max = std::max(max, (typename SendSection::point_type) (newMax+1));
710:           }
711:         }
712:         if (!rRanks->size()) {min = max = 0;}
713:         recvSection->setChart(typename RecvSection::chart_type(min, max));
714:         // TODO: deallocation
715:       };
716:       // Copy the overlap section to the related processes
717:       //   This version is for different sections, possibly with different data types
718:       // TODO: Can cache MPIMover objects (like a VecScatter)
719:       template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
720:       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) {
721:         typedef std::pair<int, typename SendSection::value_type *> allocPair;
722:         typedef typename RecvSection::point_type                   recv_point_type;
723:         const Obj<typename SendSection::atlas_type>& sendAtlas = sendSection->getAtlas();
724:         const Obj<typename RecvSection::atlas_type>& recvAtlas = recvSection->getAtlas();
725:         MPIMover<typename SendSection::value_type>   vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
726:         std::map<int, allocPair>                     sendValues;
727:         std::map<int, allocPair>                     recvValues;
728:         typename SendSection::alloc_type             sendAllocator;
729:         typename RecvSection::alloc_type             recvAllocator;

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

735:         // TODO: This should be const_iterator, but Sifter sucks
736:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
737:           const Obj<typename SendOverlap::coneSequence>&     points  = sendOverlap->cone(*r_iter);
738:           const typename SendOverlap::coneSequence::iterator pEnd    = points->end();
739:           int                                                numVals = 0;

741:           // TODO: This should be const_iterator, but Sifter sucks
742:           for(typename SendOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != pEnd; ++c_iter) {
743:             numVals += sendSection->getFiberDimension(*c_iter);
744:           }
745:           typename SendSection::value_type *v = sendAllocator.allocate(numVals);
746:           int                               k = 0;

748:           for(int i = 0; i < numVals; ++i) {sendAllocator.construct(v+i, 0);}
749:           for(typename SendOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != pEnd; ++c_iter) {
750:             const typename SendSection::value_type *vals = sendSection->restrictPoint(*c_iter);

752:             for(int i = 0; i < sendSection->getFiberDimension(*c_iter); ++i, ++k) v[k] = vals[i];
753:           }
754:           sendValues[*r_iter] = allocPair(numVals, v);
755:           vMover.send(*r_iter, numVals, v);
756:         }
757:         const Obj<typename RecvOverlap::traits::capSequence>      rRanks = recvOverlap->cap();
758:         const typename RecvOverlap::traits::capSequence::iterator rEnd   = rRanks->end();

760:         recvSection->allocatePoint();
761:         // TODO: This should be const_iterator, but Sifter sucks
762:         int maxVals = 0;
763:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
764:           const Obj<typename RecvOverlap::supportSequence>&     points  = recvOverlap->support(*r_iter);
765:           const typename RecvOverlap::supportSequence::iterator pEnd    = points->end();
766:           int                                                   numVals = 0;

768:           // TODO: This should be const_iterator, but Sifter sucks
769:           for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter) {
770:             numVals += recvSection->getFiberDimension(recv_point_type(*r_iter, s_iter.color()));
771:           }
772:           typename SendSection::value_type *v = sendAllocator.allocate(numVals);

774:           for(int i = 0; i < numVals; ++i) {sendAllocator.construct(v+i, 0);}
775:           recvValues[*r_iter] = allocPair(numVals, v);
776:           vMover.recv(*r_iter, numVals, v);
777:           maxVals = std::max(maxVals, numVals);
778:         }
779:         vMover.start();
780:         vMover.end();
781:         typename RecvSection::value_type *convertedValues = recvAllocator.allocate(maxVals);
782:         for(int i = 0; i < maxVals; ++i) {recvAllocator.construct(convertedValues+i, 0);}
783:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
784:           const Obj<typename RecvOverlap::supportSequence>&     points  = recvOverlap->support(*r_iter);
785:           const typename RecvOverlap::supportSequence::iterator pEnd    = points->end();
786:           typename SendSection::value_type                     *v       = recvValues[*r_iter].second;
787:           const int                                             numVals = recvValues[*r_iter].first;
788:           int                                                   k       = 0;

790:           for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter) {
791:             const int size = recvSection->getFiberDimension(recv_point_type(*r_iter, s_iter.color()));

793:             for(int i = 0; i < size; ++i) {convertedValues[i] = (typename RecvSection::value_type) v[k+i];}
794:             if (size) {recvSection->updatePoint(recv_point_type(*r_iter, s_iter.color()), convertedValues);}
795:             k += size;
796:           }
797:           for(int i = 0; i < numVals; ++i) {sendAllocator.destroy(v+i);}
798:         }
799:         for(int i = 0; i < maxVals; ++i) {recvAllocator.destroy(convertedValues+i);}
800:         recvAllocator.deallocate(convertedValues, maxVals);
801:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
802:           typename SendSection::value_type *v       = sendValues[*r_iter].second;
803:           const int                         numVals = sendValues[*r_iter].first;

805:           for(int i = 0; i < numVals; ++i) {sendAllocator.destroy(v+i);}
806:           sendAllocator.deallocate(v, numVals);
807:         }
808:         //recvSection->view("After copy");
809:       };
810:       // Copy the overlap section to the related processes
811:       //   This version is for sections with the same type
812:       template<typename SendOverlap, typename RecvOverlap, typename Section>
813:       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) {
814:         typedef std::pair<int, typename Section::value_type *> allocPair;
815:         const Obj<typename Section::atlas_type>& sendAtlas = sendSection->getAtlas();
816:         const Obj<typename Section::atlas_type>& recvAtlas = recvSection->getAtlas();
817:         MPIMover<typename Section::value_type>   vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
818:         std::map<int, allocPair>                 sendValues;
819:         std::map<int, allocPair>                 recvValues;
820:         typename Section::alloc_type             allocator;

822:         ///sendAtlas->view("Send Atlas in same type copy()");
823:         copy(sendOverlap, recvOverlap, sendAtlas, recvAtlas);
824:         ///recvAtlas->view("Recv Atlas after same type copy()");
825:         const Obj<typename SendOverlap::traits::baseSequence>      sRanks = sendOverlap->base();
826:         const typename SendOverlap::traits::baseSequence::iterator sEnd   = sRanks->end();

828:         // TODO: This should be const_iterator, but Sifter sucks
829:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
830:           const Obj<typename SendOverlap::coneSequence>&     points  = sendOverlap->cone(*r_iter);
831:           const typename SendOverlap::coneSequence::iterator pEnd    = points->end();
832:           int                                                numVals = 0;

834:           // TODO: This should be const_iterator, but Sifter sucks
835:           for(typename SendOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != pEnd; ++c_iter) {
836:             numVals += sendSection->getFiberDimension(*c_iter);
837:           }
838:           typename Section::value_type *v = allocator.allocate(numVals);
839:           int                           k = 0;

841:           for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
842:           for(typename SendOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != pEnd; ++c_iter) {
843:             const typename Section::value_type *vals = sendSection->restrictPoint(*c_iter);
844:             const int                           dim  = sendSection->getFiberDimension(*c_iter);

846:             for(int i = 0; i < dim; ++i, ++k) v[k] = vals[i];
847:           }
848:           sendValues[*r_iter] = allocPair(numVals, v);
849:           vMover.send(*r_iter, numVals, v);
850:         }
851:         const Obj<typename RecvOverlap::traits::capSequence>      rRanks = recvOverlap->cap();
852:         const typename RecvOverlap::traits::capSequence::iterator rEnd   = rRanks->end();

854:         recvSection->allocatePoint();
855:         ///recvSection->view("Recv Section after same type copy() allocation");
856:         // TODO: This should be const_iterator, but Sifter sucks
857:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
858:           const Obj<typename RecvOverlap::supportSequence>&     points  = recvOverlap->support(*r_iter);
859:           const typename RecvOverlap::supportSequence::iterator pEnd    = points->end();
860:           int                                                   numVals = 0;

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

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

881:           for(typename RecvOverlap::supportSequence::iterator s_iter = points->begin(); s_iter != pEnd; ++s_iter) {
882:             const int size = recvSection->getFiberDimension(s_iter.color());

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

894:           for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
895:           allocator.deallocate(v, numVals);
896:         }
897:         //recvSection->view("After copy");
898:       };
899:       // Specialize to ArrowSection
900:       template<typename SendOverlap, typename RecvOverlap>
901:       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) {
902:         typedef UniformSection<MinimalArrow<int,int>,int>      Section;
903:         typedef std::pair<int, typename Section::value_type *> allocPair;
904:         const Obj<typename Section::atlas_type>& sendAtlas = sendSection->getAtlas();
905:         const Obj<typename Section::atlas_type>& recvAtlas = recvSection->getAtlas();
906:         MPIMover<typename Section::value_type>   vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
907:         std::map<int, allocPair>                 sendValues;
908:         std::map<int, allocPair>                 recvValues;
909:         typename Section::alloc_type             allocator;

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

915:         // TODO: This should be const_iterator, but Sifter sucks
916:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
917:           const Obj<typename SendOverlap::coneSequence>&     points  = sendOverlap->cone(*r_iter);
918:           const typename SendOverlap::coneSequence::iterator pBegin  = points->begin();
919:           const typename SendOverlap::coneSequence::iterator pEnd    = points->end();
920:           int                                                numVals = 0;

922:           // TODO: This should be const_iterator, but Sifter sucks
923:           for(typename SendOverlap::coneSequence::iterator c_iter = pBegin; c_iter != pEnd; ++c_iter) {
924:             for(typename SendOverlap::coneSequence::iterator d_iter = pBegin; d_iter != pEnd; ++d_iter) {
925:               numVals += sendSection->getFiberDimension(typename Section::point_type(*c_iter, *d_iter));
926:             }
927:           }
928:           typename Section::value_type *v = allocator.allocate(numVals);
929:           int                           k = 0;

931:           for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
932:           for(typename SendOverlap::coneSequence::iterator c_iter = pBegin; c_iter != pEnd; ++c_iter) {
933:             for(typename SendOverlap::coneSequence::iterator d_iter = pBegin; d_iter != pEnd; ++d_iter) {
934:               const typename Section::point_type  arrow(*c_iter, *d_iter);
935:               const typename Section::value_type *vals = sendSection->restrictPoint(arrow);
936:               const int                           dim  = sendSection->getFiberDimension(arrow);

938:               for(int i = 0; i < dim; ++i, ++k) v[k] = vals[i];
939:             }
940:           }
941:           sendValues[*r_iter] = allocPair(numVals, v);
942:           vMover.send(*r_iter, numVals, v);
943:         }
944:         const Obj<typename RecvOverlap::traits::capSequence>      rRanks = recvOverlap->cap();
945:         const typename RecvOverlap::traits::capSequence::iterator rEnd   = rRanks->end();

947:         recvSection->allocatePoint();
948:         // TODO: This should be const_iterator, but Sifter sucks
949:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
950:           const Obj<typename RecvOverlap::supportSequence>&     points  = recvOverlap->support(*r_iter);
951:           const typename RecvOverlap::supportSequence::iterator pBegin  = points->begin();
952:           const typename RecvOverlap::supportSequence::iterator pEnd    = points->end();
953:           int                                                   numVals = 0;

955:           // TODO: This should be const_iterator, but Sifter sucks
956:           for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
957:             for(typename RecvOverlap::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter) {
958:               numVals += recvSection->getFiberDimension(typename Section::point_type(s_iter.color(), t_iter.color()));
959:             }
960:           }
961:           typename Section::value_type *v = allocator.allocate(numVals);

963:           recvValues[*r_iter] = allocPair(numVals, v);
964:           for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
965:           vMover.recv(*r_iter, numVals, v);
966:         }
967:         vMover.start();
968:         vMover.end();
969:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
970:           const Obj<typename RecvOverlap::supportSequence>&     points  = recvOverlap->support(*r_iter);
971:           const typename RecvOverlap::supportSequence::iterator pBegin  = points->begin();
972:           const typename RecvOverlap::supportSequence::iterator pEnd    = points->end();
973:           typename Section::value_type                         *v       = recvValues[*r_iter].second;
974:           const int                                             numVals = recvValues[*r_iter].first;
975:           int                                                   k       = 0;

977:           for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
978:             for(typename RecvOverlap::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter) {
979:               const typename Section::point_type arrow(s_iter.color(), t_iter.color());
980:               const int size = recvSection->getFiberDimension(arrow);

982:               if (size) {recvSection->updatePoint(arrow, &v[k]);}
983:               k += size;
984:             }
985:           }
986:           for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
987:           allocator.deallocate(v, numVals);
988:         }
989:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
990:           typename Section::value_type *v       = sendValues[*r_iter].second;
991:           const int                     numVals = sendValues[*r_iter].first;

993:           for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
994:           allocator.deallocate(v, numVals);
995:         }
996:         //recvSection->view("After copy");
997:       };
998:       // Specialize to a ConstantSection
999: #if 0
1000:       template<typename SendOverlap, typename RecvOverlap, typename Value>
1001:       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) {
1002:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1003:       };
1004:       template<typename SendOverlap, typename RecvOverlap, typename Value>
1005:       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) {
1006:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1007:       };
1008: #else
1009:       template<typename SendOverlap, typename RecvOverlap, typename Value>
1010:       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) {
1011:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1012:       };
1013:       template<typename SendOverlap, typename RecvOverlap, typename Value>
1014:       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) {
1015:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1016:       };
1017: #endif
1018:       // Specialize to an IConstantSection
1019:       template<typename SendOverlap, typename RecvOverlap, typename Value>
1020:       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) {
1021:         // Why doesn't this work?
1022:         //   This supposed to be a copy, BUT filtered through the sendOverlap
1023:         //   However, an IConstant section does not allow filtration of its
1024:         //   chart. Therefore, you end up with either
1025:         //
1026:         //   a) Too many items in the chart, copied from the remote sendSection
1027:         //   b) A chart mapped to the local numbering, which we do not want
1028:         copyIConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1029:       };
1030:       // Specialize to an BaseSection/ConstantSection pair
1031: #if 0
1032:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
1033:       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) {
1034:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1035:       };
1036: #else
1037:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
1038:       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) {
1039:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1040:       };
1041: #endif
1042:       // Specialize to an BaseSectionV/ConstantSection pair
1043: #if 0
1044:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
1045:       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) {
1046:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1047:       };
1048: #else
1049:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
1050:       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) {
1051:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1052:       };
1053: #endif
1054:       // Specialize to an LabelBaseSection/ConstantSection pair
1055: #if 0
1056:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_, typename Label_>
1057:       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) {
1058:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1059:       };
1060: #else
1061:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_, typename Label_>
1062:       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) {
1063:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1064:       };
1065: #endif
1066:       // Specialize to a ConstantSection for ArrowSection
1067:       template<typename SendOverlap, typename RecvOverlap, typename Value>
1068:       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) {
1069:         copyConstantArrow(sendOverlap, recvOverlap, sendSection, recvSection);
1070:       };
1071:     };
1072:     class BinaryFusion {
1073:     public:
1074:       template<typename OverlapSection, typename RecvOverlap, typename Section, typename Function>
1075:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section, Function binaryOp) {
1076:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
1077:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

1079:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1080:           // TODO: This must become a more general iterator over colors
1081:           const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1082:           // Just taking the first value
1083:           const typename Section::point_type&        localPoint    = *p_iter;
1084:           const typename OverlapSection::point_type& remotePoint   = points->begin().color();
1085:           const typename OverlapSection::value_type *overlapValues = overlapSection->restrictPoint(remotePoint);
1086:           const typename Section::value_type        *localValues   = section->restrictPoint(localPoint);
1087:           const int                                  dim           = section->getFiberDimension(localPoint);
1088:           // TODO: optimize allocation
1089:           typename Section::value_type              *values        = new typename Section::value_type[dim];

1091:           for(int d = 0; d < dim; ++d) {
1092:             values[d] = binaryOp(overlapValues[d], localValues[d]);
1093:           }
1094:           section->updatePoint(localPoint, values);
1095:           delete [] values;
1096:         }
1097:       };
1098:     };
1099:     class ReplacementBinaryFusion {
1100:     public:
1101:       template<typename OverlapSection, typename RecvOverlap, typename Section>
1102:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
1103:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
1104:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

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

1113:           section->update(localPoint, overlapSection->restrictPoint(remotePoint));
1114:         }
1115:       };
1116:     };
1117:     class AdditiveBinaryFusion {
1118:     public:
1119:       template<typename OverlapSection, typename RecvOverlap, typename Section>
1120:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
1121:         typedef typename Section::point_type        point_type;
1122:         typedef typename OverlapSection::point_type overlap_point_type;
1123:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
1124:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

1126:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1127:           // TODO: This must become a more general iterator over colors
1128:           const typename Section::point_type&                localPoint = *p_iter;
1129:           const Obj<typename RecvOverlap::coneSequence>&     points     = recvOverlap->cone(*p_iter);
1130:           const typename RecvOverlap::coneSequence::iterator cEnd       = points->end();

1132:           for(typename RecvOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != cEnd; ++c_iter) {
1133:             const int         rank        = *c_iter;
1134:             const point_type& remotePoint = c_iter.color();

1136:             section->updateAddPoint(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
1137:           }
1138:         }
1139:       };
1140:     };
1141:     class InsertionBinaryFusion {
1142:     public:
1143:       // Insert the overlapSection values into section along recvOverlap
1144:       template<typename OverlapSection, typename RecvOverlap, typename Section>
1145:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
1146:         typedef typename Section::point_type        point_type;
1147:         typedef typename OverlapSection::point_type overlap_point_type;
1148:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
1149:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

1151:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1152:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
1153:           const point_type&                              localPoint  = *p_iter;
1154:           const int                                      rank        = *points->begin();
1155:           const point_type&                              remotePoint = points->begin().color();

1157:           if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1158:             if (!section->hasPoint(localPoint)) {
1159:               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;
1160:             }
1161:             section->setFiberDimension(localPoint, overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)));
1162:           }
1163:         }
1164:         if (rPoints->size()) {section->allocatePoint();}
1165:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1166:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
1167:           const point_type&                              localPoint  = *p_iter;
1168:           const int                                      rank        = *points->begin();
1169:           const point_type&                              remotePoint = points->begin().color();

1171:           if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1172:             section->updatePoint(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
1173:           }
1174:         }
1175:       };
1176:       // Specialize to ArrowSection
1177:       template<typename OverlapSection, typename RecvOverlap>
1178:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<UniformSection<MinimalArrow<int,int>,int> >& section) {
1179:         typedef UniformSection<MinimalArrow<int,int>,int> Section;
1180:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
1181:         const typename RecvOverlap::traits::baseSequence::iterator rBegin  = rPoints->begin();
1182:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

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

1189:           for(typename RecvOverlap::traits::baseSequence::iterator q_iter = rBegin; q_iter != rEnd; ++q_iter) {
1190:             const Obj<typename RecvOverlap::coneSequence>& targets      = recvOverlap->cone(*q_iter);
1191:             const typename RecvOverlap::target_type        localTarget  = *q_iter;
1192:             const typename RecvOverlap::target_type        remoteTarget = targets->begin().color();
1193:             const typename Section::point_type             localPoint(localSource, localTarget);
1194:             const typename OverlapSection::point_type      remotePoint(remoteSource, remoteTarget);

1196:             if (overlapSection->hasPoint(remotePoint)) {section->setFiberDimension(localPoint, overlapSection->getFiberDimension(remotePoint));}
1197:           }
1198:         }
1199:         if (rPoints->size()) {section->allocatePoint();}
1200:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rBegin; p_iter != rEnd; ++p_iter) {
1201:           const Obj<typename RecvOverlap::coneSequence>& sources      = recvOverlap->cone(*p_iter);
1202:           const typename RecvOverlap::target_type        localSource  = *p_iter;
1203:           const typename RecvOverlap::target_type        remoteSource = sources->begin().color();

1205:           for(typename RecvOverlap::traits::baseSequence::iterator q_iter = rBegin; q_iter != rEnd; ++q_iter) {
1206:             const Obj<typename RecvOverlap::coneSequence>& targets      = recvOverlap->cone(*q_iter);
1207:             const typename RecvOverlap::target_type        localTarget  = *q_iter;
1208:             const typename RecvOverlap::target_type        remoteTarget = targets->begin().color();
1209:             const typename Section::point_type             localPoint(localSource, localTarget);
1210:             const typename OverlapSection::point_type      remotePoint(remoteSource, remoteTarget);
1211: 
1212:             if (overlapSection->hasPoint(remotePoint)) {section->updatePoint(localPoint, overlapSection->restrictPoint(remotePoint));}
1213:           }
1214:         }
1215:       };
1216:       // Specialize to the Sieve
1217:       template<typename OverlapSection, typename RecvOverlap, typename Renumbering, typename Point>
1218:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, Renumbering& renumbering, const Obj<Sieve<Point,Point,int> >& sieve) {
1219:         typedef typename OverlapSection::point_type overlap_point_type;
1220:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
1221:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

1223:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1224:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
1225:           const Point&                                   localPoint  = *p_iter;
1226:           const int                                      rank        = *points->begin();
1227:           const Point&                                   remotePoint = points->begin().color();
1228:           const int                                      size        = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1229:           const typename OverlapSection::value_type     *values      = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1230:           int                                            c           = 0;

1232:           sieve->clearCone(localPoint);
1233:           for(int i = 0; i < size; ++i, ++c) {sieve->addCone(renumbering[values[i]], localPoint, c);}
1234:         }
1235:       };
1236:       // Specialize to the ISieve
1237:       template<typename OverlapSection, typename RecvOverlap, typename Renumbering, typename Point>
1238:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, Renumbering& renumbering, const Obj<IFSieve<Point> >& sieve) {
1239:         typedef typename OverlapSection::point_type overlap_point_type;
1240:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
1241:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();
1242:         int                                                        maxSize = 0;

1244:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1245:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
1246:           const Point&                                   localPoint  = *p_iter;
1247:           const int                                      rank        = *points->begin();
1248:           const Point&                                   remotePoint = points->begin().color();
1249:           const int                                      size        = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1250:           const typename OverlapSection::value_type     *values      = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));

1252:           sieve->setConeSize(localPoint, size);
1253:           ///for(int i = 0; i < size; ++i) {sieve->addSupportSize(renumbering[values[i]], 1);}
1254:           for(int i = 0; i < size; ++i) {sieve->addSupportSize(renumbering[values[i].first], 1);}
1255:           maxSize = std::max(maxSize, size);
1256:         }
1257:         sieve->allocate();
1258:         ///typename OverlapSection::value_type *localValues = new typename OverlapSection::value_type[maxSize];
1259:         typename OverlapSection::value_type::first_type  *localValues      = new typename OverlapSection::value_type::first_type[maxSize];
1260:         typename OverlapSection::value_type::second_type *localOrientation = new typename OverlapSection::value_type::second_type[maxSize];

1262:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1263:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
1264:           const Point&                                   localPoint  = *p_iter;
1265:           const int                                      rank        = *points->begin();
1266:           const Point&                                   remotePoint = points->begin().color();
1267:           const int                                      size        = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1268:           const typename OverlapSection::value_type     *values      = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));

1270:           ///for(int i = 0; i < size; ++i) {localValues[i] = renumbering[values[i]];}
1271:           for(int i = 0; i < size; ++i) {
1272:             localValues[i]      = renumbering[values[i].first];
1273:             localOrientation[i] = values[i].second;
1274:           }
1275:           sieve->setCone(localValues, localPoint);
1276:           sieve->setConeOrientation(localOrientation, localPoint);
1277:         }
1278:         delete [] localValues;
1279:         delete [] localOrientation;
1280:       };
1281:       template<typename OverlapSection, typename RecvOverlap, typename Point>
1282:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<IFSieve<Point> >& sieve) {
1283:         typedef typename OverlapSection::point_type overlap_point_type;
1284:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
1285:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();
1286:         int                                                        maxSize = 0;

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

1296:           sieve->setConeSize(localPoint, size);
1297:           for(int i = 0; i < size; ++i) {sieve->addSupportSize(values[i].first, 1);}
1298:           maxSize = std::max(maxSize, size);
1299:         }
1300:         sieve->allocate();
1301:         typename OverlapSection::value_type::first_type  *localValues      = new typename OverlapSection::value_type::first_type[maxSize];
1302:         typename OverlapSection::value_type::second_type *localOrientation = new typename OverlapSection::value_type::second_type[maxSize];

1304:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1305:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
1306:           const Point&                                   localPoint  = *p_iter;
1307:           const int                                      rank        = *points->begin();
1308:           const Point&                                   remotePoint = points->begin().color();
1309:           const int                                      size        = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1310:           const typename OverlapSection::value_type     *values      = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));

1312:           for(int i = 0; i < size; ++i) {
1313:             localValues[i]      = values[i].first;
1314:             localOrientation[i] = values[i].second;
1315:           }
1316:           sieve->setCone(localValues, localPoint);
1317:           sieve->setConeOrientation(localOrientation, localPoint);
1318:         }
1319:         delete [] localValues;
1320:         delete [] localOrientation;
1321:       };
1322:       // Generic
1323:       template<typename OverlapSection, typename RecvOverlap, typename Section, typename Bundle>
1324:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section, const Obj<Bundle>& bundle) {
1325:         typedef typename OverlapSection::point_type overlap_point_type;
1326:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
1327:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

1329:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1330:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
1331:           const typename Section::point_type&            localPoint  = *p_iter;
1332:           const int                                      rank        = *points->begin();
1333:           const typename OverlapSection::point_type&     remotePoint = points->begin().color();

1335:           section->setFiberDimension(localPoint, overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)));
1336:         }
1337:         bundle->allocate(section);
1338:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1339:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
1340:           const typename Section::point_type&            localPoint  = *p_iter;
1341:           const int                                      rank        = *points->begin();
1342:           const typename OverlapSection::point_type&     remotePoint = points->begin().color();

1344:           section->update(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
1345:         }
1346:       };
1347:     };
1348:   }
1349: }

1351: #endif