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:         int MPI_Type_free(&this->_datatype);CHKERRXX(ierr);
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], ierr;
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);CHKERRXX(ierr);
215:         MPI_Type_commit(&newtype);CHKERRXX(ierr);
216:         this->_createdType = true;
217:         return newtype;
218:       } else if (sizeof(value_type) == 32) {
219:         int          blen[2], ierr;
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);CHKERRXX(ierr);
225:         MPI_Type_commit(&newtype);CHKERRXX(ierr);
226:         this->_createdType = true;
227:         return newtype;
228:       }
229:       ostringstream msg;

231:       msg << "Cannot determine MPI type for value type with size " << sizeof(value_type);
232:       throw PETSc::Exception(msg.str().c_str());
233:     };
234:     int getNewTag() const {
235:       static int tagKeyval = MPI_KEYVAL_INVALID;
236:       int *tagvalp = NULL, *maxval, flg, ierr;

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

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

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

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

302:         int MPI_Start(&request);CHKERRXX(ierr);
303:       }
304:     };
305:     void end() {
306:       MPI_Status status;

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

311:         int MPI_Wait(&request, &status);CHKERRXX(ierr);
312:       }
313:       for(typename requests_type::const_iterator r_iter = this->_requests.begin(); r_iter != this->_requests.end(); ++r_iter) {
314:         MPI_Request request = *r_iter;

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

387:       // Change to Allgather() for the correct binning algorithm
388:       MPI_Gather(&size, 1, MPI_INT, sizes, 1, MPI_INT, 0, comm);
389:       if (commRank == 0) {
390:         offsets[0] = 0;
391:         for(int p = 1; p <= commSize; p++) {
392:           offsets[p] = offsets[p-1] + sizes[p-1];
393:         }
394:         numRemotePoints = offsets[commSize];
395:         remotePoints    = pointPair_allocator.allocate(numRemotePoints/2);
396:         for(int i = 0; i < numRemotePoints/2; ++i) {pointPair_allocator.construct(remotePoints+i, pointPair());}
397:       }
398:       MPI_Gatherv(sendBuf, size, MPI_INT, remotePoints, sizes, offsets, MPI_INT, 0, comm);
399:       for(size_t i = 0; i < points.size(); ++i) {point_allocator.destroy(sendBuf+i);}
400:       point_allocator.deallocate(sendBuf, points.size());
401:       std::map<int, std::map<int, std::set<pointTriple> > > overlapInfo; // Maps (p,q) to their set of overlap points

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

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

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

481:           if (debug) {std::cout << "["<<commRank<<"]Matched up remote point " << remotePoint << "("<<point<<") to local " << renumbering[point] << std::endl;}
482:           sendOverlap->addArrow(renumbering[point], rank, remotePoint);
483:           recvOverlap->addArrow(rank, renumbering[point], remotePoint);
484:         }
485:       }

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

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

521:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
522:           const Obj<typename SendOverlap::coneSequence>&     points = sendOverlap->cone(*r_iter);
523:           const int                                          pSize  = points->size();
524:           const typename SendOverlap::coneSequence::iterator pEnd   = points->end();
525:           char                                              *v      = sendAllocator.allocate(points->size());
526:           int                                                k      = 0;

528:           for(int i = 0; i < pSize; ++i) {sendAllocator.construct(v+i, 0);}
529:           for(typename SendOverlap::coneSequence::iterator p_iter = points->begin(); p_iter != pEnd; ++p_iter, ++k) {
530:             v[k] = (char) sendSection->hasPoint(*p_iter);
531:           }
532:           sendPoints[*r_iter] = v;
533:           pMover.send(*r_iter, pSize, sendPoints[*r_iter]);
534:           vMover.send(*r_iter, 2, sValues);
535:         }
536:         const Obj<typename RecvOverlap::traits::capSequence>      rRanks  = recvOverlap->cap();
537:         const typename RecvOverlap::traits::capSequence::iterator rEnd    = rRanks->end();
538:         const typename RecvSection::value_type                   *rValues = recvSection->restrictSpace();

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

545:           for(int i = 0; i < pSize; ++i) {recvAllocator.construct(v+i, 0);}
546:           recvPoints[*r_iter] = v;
547:           pMover.recv(*r_iter, pSize, recvPoints[*r_iter]);
548:           vMover.recv(*r_iter, 2, rValues);
549:         }
550:         pMover.start();
551:         pMover.end();
552:         vMover.start();
553:         vMover.end();
554:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
555:           const Obj<typename RecvOverlap::traits::supportSequence>&     points = recvOverlap->support(*r_iter);
556:           const typename RecvOverlap::traits::supportSequence::iterator pEnd   = points->end();
557:           const char                                                   *v      = recvPoints[*r_iter];
558:           int                                                           k      = 0;

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

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

585:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
586:           const Obj<typename SendOverlap::coneSequence>&     points = sendOverlap->cone(*r_iter);
587:           const int                                          pSize  = points->size();
588:           const typename SendOverlap::coneSequence::iterator pBegin = points->begin();
589:           const typename SendOverlap::coneSequence::iterator pEnd   = points->end();
590:           char                                              *v      = sendAllocator.allocate(pSize*pSize);
591:           int                                                k      = 0;

593:           for(size_t i = 0; i < pSize*pSize; ++i) {sendAllocator.construct(v+i, 0);}
594:           for(typename SendOverlap::coneSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
595:             for(typename SendOverlap::coneSequence::iterator q_iter = pBegin; q_iter != pEnd; ++q_iter, ++k) {
596:               v[k] = (char) sendSection->hasPoint(typename SendSection::point_type(*p_iter,*q_iter));
597:             }
598:           }
599:           sendPoints[*r_iter] = v;
600:           pMover.send(*r_iter, pSize*pSize, sendPoints[*r_iter]);
601:           vMover.send(*r_iter, 2, sValues);
602:         }
603:         const Obj<typename RecvOverlap::traits::capSequence>      rRanks  = recvOverlap->cap();
604:         const typename RecvOverlap::traits::capSequence::iterator rEnd    = rRanks->end();
605:         const typename RecvSection::value_type                   *rValues = recvSection->restrictSpace();

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

612:           for(size_t i = 0; i < pSize*pSize; ++i) {recvAllocator.construct(v+i, 0);}
613:           recvPoints[*r_iter] = v;
614:           pMover.recv(*r_iter, pSize*pSize, recvPoints[*r_iter]);
615:           vMover.recv(*r_iter, 2, rValues);
616:         }
617:         pMover.start();
618:         pMover.end();
619:         vMover.start();
620:         vMover.end();
621:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
622:           const Obj<typename RecvOverlap::traits::supportSequence>&     points = recvOverlap->support(*r_iter);
623:           const typename RecvOverlap::traits::supportSequence::iterator pBegin = points->begin();
624:           const typename RecvOverlap::traits::supportSequence::iterator pEnd   = points->end();
625:           const char                                                   *v      = recvPoints[*r_iter];
626:           int                                                           k      = 0;

628:           for(typename RecvOverlap::traits::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
629:             for(typename RecvOverlap::traits::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter, ++k) {
630:               if (v[k]) {recvSection->addPoint(typename RecvSection::point_type(s_iter.color(),t_iter.color()));}
631:             }
632:           }
633:         }
634:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
635:           sendAllocator.deallocate(sendPoints[*r_iter], sendOverlap->cone(*r_iter)->size()*sendOverlap->cone(*r_iter)->size());
636:         }
637:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
638:           recvAllocator.deallocate(recvPoints[*r_iter], recvOverlap->support(*r_iter)->size()*recvOverlap->support(*r_iter)->size());
639:         }
640:       };
641:       // Copy the overlap section to the related processes
642:       //   This version is for IConstant sections, meaning the same, single value over all points
643:       template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
644:       static void copyIConstant(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
645:         MPIMover<typename SendSection::point_type> pMover(sendSection->comm(), sendSection->debug());
646:         MPIMover<typename SendSection::value_type> vMover(sendSection->comm(), sendSection->debug());
647:         std::map<int, typename SendSection::point_type *> sendPoints;
648:         std::map<int, typename SendSection::point_type *> recvPoints;
649:         typename SendSection::alloc_type::template rebind<typename SendSection::point_type>::other sendAllocator;
650:         typename RecvSection::alloc_type::template rebind<typename SendSection::point_type>::other recvAllocator;

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

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

659:           for(size_t i = 0; i < 2; ++i) {sendAllocator.construct(v+i, 0);}
660:           v[0] = sendSection->getChart().min();
661:           v[1] = sendSection->getChart().max();
662:           sendPoints[*r_iter] = v;
663:           pMover.send(*r_iter, 2, sendPoints[*r_iter]);
664:           vMover.send(*r_iter, 2, sValues);
665:           ///std::cout << "["<<sendOverlap->commRank()<<"]Sending chart (" << v[0] << ", " << v[1] << ") with values (" << sValues[0] << ", " << sValues[1] << ") to process " << *r_iter << std::endl;
666:         }
667:         const Obj<typename RecvOverlap::traits::capSequence>      rRanks  = recvOverlap->cap();
668:         const typename RecvOverlap::traits::capSequence::iterator rEnd    = rRanks->end();
669:         const typename RecvSection::value_type                   *rValues = recvSection->restrictSpace();

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

674:           for(size_t i = 0; i < 2; ++i) {recvAllocator.construct(v+i, 0);}
675:           recvPoints[*r_iter] = v;
676:           pMover.recv(*r_iter, 2, recvPoints[*r_iter]);
677:           vMover.recv(*r_iter, 2, rValues);
678:         }
679:         pMover.start();
680:         pMover.end();
681:         vMover.start();
682:         vMover.end();

684:         typename SendSection::point_type min = -1;
685:         typename SendSection::point_type max = -1;

687:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
688:           const typename RecvSection::point_type *v = recvPoints[*r_iter];
689:           typename SendSection::point_type        newMin = v[0];
690:           typename SendSection::point_type        newMax = v[1]-1;
691:           //int                                     pSize  = 0;

693:           ///std::cout << "["<<recvOverlap->commRank()<<"]Received chart (" << v[0] << ", " << v[1] << ") from process " << *r_iter << std::endl;
694: #if 0
695:           // Translate to local numbering
696:           if (recvOverlap->support(*r_iter)->size()) {
697:             while(!pSize) {
698:               const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter, newMin);
699:               pSize = points->size();
700:               if (pSize) {
701:                 newMin = *points->begin();
702:               } else {
703:                 newMin++;
704:               }
705:             }
706:             pSize  = 0;
707:             while(!pSize) {
708:               const Obj<typename RecvOverlap::supportSequence>& points = recvOverlap->support(*r_iter, newMax);
709:               pSize = points->size();
710:               if (pSize) {
711:                 newMax = *points->begin();
712:               } else {
713:                 newMax--;
714:               }
715:             }
716:           }
717:           std::cout << "["<<recvOverlap->commRank()<<"]Translated to chart (" << newMin << ", " << newMax+1 << ") from process " << *r_iter << std::endl;
718: #endif
719:           // Update chart
720:           if (min < 0) {
721:             min = newMin;
722:             max = newMax+1;
723:           } else {
724:             min = std::min(min, newMin);
725:             max = std::max(max, (typename SendSection::point_type) (newMax+1));
726:           }
727:         }
728:         if (!rRanks->size()) {min = max = 0;}
729:         recvSection->setChart(typename RecvSection::chart_type(min, max));
730:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
731:           sendAllocator.deallocate(sendPoints[*r_iter], 2);
732:         }
733:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
734:           recvAllocator.deallocate(recvPoints[*r_iter], 2);
735:         }
736:       };
737:       // Copy the overlap section to the related processes
738:       //   This version is for different sections, possibly with different data types
739:       // TODO: Can cache MPIMover objects (like a VecScatter)
740:       template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
741:       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) {
742:         typedef std::pair<int, typename SendSection::value_type *> allocPair;
743:         typedef typename RecvSection::point_type                   recv_point_type;
744:         const Obj<typename SendSection::atlas_type>& sendAtlas = sendSection->getAtlas();
745:         const Obj<typename RecvSection::atlas_type>& recvAtlas = recvSection->getAtlas();
746:         MPIMover<typename SendSection::value_type>   vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
747:         std::map<int, allocPair>                     sendValues;
748:         std::map<int, allocPair>                     recvValues;
749:         typename SendSection::alloc_type             sendAllocator;
750:         typename RecvSection::alloc_type             recvAllocator;

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

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

762:           // TODO: This should be const_iterator, but Sifter sucks
763:           for(typename SendOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != pEnd; ++c_iter) {
764:             numVals += sendSection->getFiberDimension(*c_iter);
765:           }
766:           typename SendSection::value_type *v = sendAllocator.allocate(numVals);
767:           int                               k = 0;

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

773:             for(int i = 0; i < sendSection->getFiberDimension(*c_iter); ++i, ++k) v[k] = vals[i];
774:           }
775:           sendValues[*r_iter] = allocPair(numVals, v);
776:           vMover.send(*r_iter, numVals, v);
777:         }
778:         const Obj<typename RecvOverlap::traits::capSequence>      rRanks = recvOverlap->cap();
779:         const typename RecvOverlap::traits::capSequence::iterator rEnd   = rRanks->end();

781:         recvSection->allocatePoint();
782:         // TODO: This should be const_iterator, but Sifter sucks
783:         int maxVals = 0;
784:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
785:           const Obj<typename RecvOverlap::supportSequence>&     points  = recvOverlap->support(*r_iter);
786:           const typename RecvOverlap::supportSequence::iterator pEnd    = points->end();
787:           int                                                   numVals = 0;

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

795:           for(int i = 0; i < numVals; ++i) {sendAllocator.construct(v+i, 0);}
796:           recvValues[*r_iter] = allocPair(numVals, v);
797:           vMover.recv(*r_iter, numVals, v);
798:           maxVals = std::max(maxVals, numVals);
799:         }
800:         vMover.start();
801:         vMover.end();
802:         typename RecvSection::value_type *convertedValues = recvAllocator.allocate(maxVals);
803:         for(int i = 0; i < maxVals; ++i) {recvAllocator.construct(convertedValues+i, 0);}
804:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
805:           const Obj<typename RecvOverlap::supportSequence>&     points  = recvOverlap->support(*r_iter);
806:           const typename RecvOverlap::supportSequence::iterator pEnd    = points->end();
807:           typename SendSection::value_type                     *v       = recvValues[*r_iter].second;
808:           const int                                             numVals = recvValues[*r_iter].first;
809:           int                                                   k       = 0;

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

814:             for(int i = 0; i < size; ++i) {convertedValues[i] = (typename RecvSection::value_type) v[k+i];}
815:             if (size) {recvSection->updatePoint(recv_point_type(*r_iter, s_iter.color()), convertedValues);}
816:             k += size;
817:           }
818:           for(int i = 0; i < numVals; ++i) {sendAllocator.destroy(v+i);}
819:         }
820:         for(int i = 0; i < maxVals; ++i) {recvAllocator.destroy(convertedValues+i);}
821:         recvAllocator.deallocate(convertedValues, maxVals);
822:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
823:           typename SendSection::value_type *v       = sendValues[*r_iter].second;
824:           const int                         numVals = sendValues[*r_iter].first;

826:           for(int i = 0; i < numVals; ++i) {sendAllocator.destroy(v+i);}
827:           sendAllocator.deallocate(v, numVals);
828:         }
829:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
830:           typename SendSection::value_type *v       = recvValues[*r_iter].second;
831:           const int                         numVals = recvValues[*r_iter].first;

833:           for(int i = 0; i < numVals; ++i) {sendAllocator.destroy(v+i);}
834:           sendAllocator.deallocate(v, numVals);
835:         }
836:         //recvSection->view("After copy");
837:       };
838:       // Copy the overlap section to the related processes
839:       //   This version is for sections with the same type
840:       template<typename SendOverlap, typename RecvOverlap, typename Section>
841:       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) {
842:         typedef std::pair<int, typename Section::value_type *> allocPair;
843:         const Obj<typename Section::atlas_type>& sendAtlas = sendSection->getAtlas();
844:         const Obj<typename Section::atlas_type>& recvAtlas = recvSection->getAtlas();
845:         MPIMover<typename Section::value_type>   vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
846:         std::map<int, allocPair>                 sendValues;
847:         std::map<int, allocPair>                 recvValues;
848:         typename Section::alloc_type             allocator;

850:         ///sendAtlas->view("Send Atlas in same type copy()");
851:         copy(sendOverlap, recvOverlap, sendAtlas, recvAtlas);
852:         ///recvAtlas->view("Recv Atlas after same type copy()");
853:         const Obj<typename SendOverlap::traits::baseSequence>      sRanks = sendOverlap->base();
854:         const typename SendOverlap::traits::baseSequence::iterator sEnd   = sRanks->end();

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

862:           // TODO: This should be const_iterator, but Sifter sucks
863:           for(typename SendOverlap::coneSequence::iterator c_iter = points->begin(); c_iter != pEnd; ++c_iter) {
864:             numVals += sendSection->getFiberDimension(*c_iter);
865:           }
866:           typename Section::value_type *v = allocator.allocate(numVals);
867:           int                           k = 0;

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

874:             for(int i = 0; i < dim; ++i, ++k) v[k] = vals[i];
875:           }
876:           sendValues[*r_iter] = allocPair(numVals, v);
877:           vMover.send(*r_iter, numVals, v);
878:         }
879:         const Obj<typename RecvOverlap::traits::capSequence>      rRanks = recvOverlap->cap();
880:         const typename RecvOverlap::traits::capSequence::iterator rEnd   = rRanks->end();

882:         recvSection->allocatePoint();
883:         ///recvSection->view("Recv Section after same type copy() allocation");
884:         // TODO: This should be const_iterator, but Sifter sucks
885:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
886:           const Obj<typename RecvOverlap::supportSequence>&     points  = recvOverlap->support(*r_iter);
887:           const typename RecvOverlap::supportSequence::iterator pEnd    = points->end();
888:           int                                                   numVals = 0;

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

896:           recvValues[*r_iter] = allocPair(numVals, v);
897:           for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
898:           vMover.recv(*r_iter, numVals, v);
899:         }
900:         vMover.start();
901:         vMover.end();
902:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
903:           const Obj<typename RecvOverlap::supportSequence>&     points  = recvOverlap->support(*r_iter);
904:           const typename RecvOverlap::supportSequence::iterator pEnd    = points->end();
905:           typename Section::value_type                         *v       = recvValues[*r_iter].second;
906:           const int                                             numVals = recvValues[*r_iter].first;
907:           int                                                   k       = 0;

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

912:             if (size) {recvSection->updatePoint(s_iter.color(), &v[k]);}
913:             k += size;
914:           }
915:           for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
916:           allocator.deallocate(v, numVals);
917:         }
918:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
919:           typename Section::value_type *v       = sendValues[*r_iter].second;
920:           const int                     numVals = sendValues[*r_iter].first;

922:           for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
923:           allocator.deallocate(v, numVals);
924:         }
925:         //recvSection->view("After copy");
926:       };
927:       // Specialize to ArrowSection
928:       template<typename SendOverlap, typename RecvOverlap>
929:       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) {
930:         typedef UniformSection<MinimalArrow<int,int>,int>      Section;
931:         typedef std::pair<int, typename Section::value_type *> allocPair;
932:         const Obj<typename Section::atlas_type>& sendAtlas = sendSection->getAtlas();
933:         const Obj<typename Section::atlas_type>& recvAtlas = recvSection->getAtlas();
934:         MPIMover<typename Section::value_type>   vMover(sendSection->comm(), datatype, MPI_UNDEFINED, sendSection->debug());
935:         std::map<int, allocPair>                 sendValues;
936:         std::map<int, allocPair>                 recvValues;
937:         typename Section::alloc_type             allocator;

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

943:         // TODO: This should be const_iterator, but Sifter sucks
944:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
945:           const Obj<typename SendOverlap::coneSequence>&     points  = sendOverlap->cone(*r_iter);
946:           const typename SendOverlap::coneSequence::iterator pBegin  = points->begin();
947:           const typename SendOverlap::coneSequence::iterator pEnd    = points->end();
948:           int                                                numVals = 0;

950:           // TODO: This should be const_iterator, but Sifter sucks
951:           for(typename SendOverlap::coneSequence::iterator c_iter = pBegin; c_iter != pEnd; ++c_iter) {
952:             for(typename SendOverlap::coneSequence::iterator d_iter = pBegin; d_iter != pEnd; ++d_iter) {
953:               numVals += sendSection->getFiberDimension(typename Section::point_type(*c_iter, *d_iter));
954:             }
955:           }
956:           typename Section::value_type *v = allocator.allocate(numVals);
957:           int                           k = 0;

959:           for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
960:           for(typename SendOverlap::coneSequence::iterator c_iter = pBegin; c_iter != pEnd; ++c_iter) {
961:             for(typename SendOverlap::coneSequence::iterator d_iter = pBegin; d_iter != pEnd; ++d_iter) {
962:               const typename Section::point_type  arrow(*c_iter, *d_iter);
963:               const typename Section::value_type *vals = sendSection->restrictPoint(arrow);
964:               const int                           dim  = sendSection->getFiberDimension(arrow);

966:               for(int i = 0; i < dim; ++i, ++k) v[k] = vals[i];
967:             }
968:           }
969:           sendValues[*r_iter] = allocPair(numVals, v);
970:           vMover.send(*r_iter, numVals, v);
971:         }
972:         const Obj<typename RecvOverlap::traits::capSequence>      rRanks = recvOverlap->cap();
973:         const typename RecvOverlap::traits::capSequence::iterator rEnd   = rRanks->end();

975:         recvSection->allocatePoint();
976:         // TODO: This should be const_iterator, but Sifter sucks
977:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
978:           const Obj<typename RecvOverlap::supportSequence>&     points  = recvOverlap->support(*r_iter);
979:           const typename RecvOverlap::supportSequence::iterator pBegin  = points->begin();
980:           const typename RecvOverlap::supportSequence::iterator pEnd    = points->end();
981:           int                                                   numVals = 0;

983:           // TODO: This should be const_iterator, but Sifter sucks
984:           for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
985:             for(typename RecvOverlap::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter) {
986:               numVals += recvSection->getFiberDimension(typename Section::point_type(s_iter.color(), t_iter.color()));
987:             }
988:           }
989:           typename Section::value_type *v = allocator.allocate(numVals);

991:           recvValues[*r_iter] = allocPair(numVals, v);
992:           for(int i = 0; i < numVals; ++i) {allocator.construct(v+i, 0);}
993:           vMover.recv(*r_iter, numVals, v);
994:         }
995:         vMover.start();
996:         vMover.end();
997:         for(typename RecvOverlap::traits::capSequence::iterator r_iter = rRanks->begin(); r_iter != rEnd; ++r_iter) {
998:           const Obj<typename RecvOverlap::supportSequence>&     points  = recvOverlap->support(*r_iter);
999:           const typename RecvOverlap::supportSequence::iterator pBegin  = points->begin();
1000:           const typename RecvOverlap::supportSequence::iterator pEnd    = points->end();
1001:           typename Section::value_type                         *v       = recvValues[*r_iter].second;
1002:           const int                                             numVals = recvValues[*r_iter].first;
1003:           int                                                   k       = 0;

1005:           for(typename RecvOverlap::supportSequence::iterator s_iter = pBegin; s_iter != pEnd; ++s_iter) {
1006:             for(typename RecvOverlap::supportSequence::iterator t_iter = pBegin; t_iter != pEnd; ++t_iter) {
1007:               const typename Section::point_type arrow(s_iter.color(), t_iter.color());
1008:               const int size = recvSection->getFiberDimension(arrow);

1010:               if (size) {recvSection->updatePoint(arrow, &v[k]);}
1011:               k += size;
1012:             }
1013:           }
1014:           for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
1015:           allocator.deallocate(v, numVals);
1016:         }
1017:         for(typename SendOverlap::traits::baseSequence::iterator r_iter = sRanks->begin(); r_iter != sEnd; ++r_iter) {
1018:           typename Section::value_type *v       = sendValues[*r_iter].second;
1019:           const int                     numVals = sendValues[*r_iter].first;

1021:           for(int i = 0; i < numVals; ++i) {allocator.destroy(v+i);}
1022:           allocator.deallocate(v, numVals);
1023:         }
1024:         //recvSection->view("After copy");
1025:       };
1026:       // Specialize to a ConstantSection
1027: #if 0
1028:       template<typename SendOverlap, typename RecvOverlap, typename Value>
1029:       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) {
1030:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1031:       };
1032:       template<typename SendOverlap, typename RecvOverlap, typename Value>
1033:       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) {
1034:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1035:       };
1036: #else
1037:       template<typename SendOverlap, typename RecvOverlap, typename Value>
1038:       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) {
1039:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1040:       };
1041:       template<typename SendOverlap, typename RecvOverlap, typename Value>
1042:       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) {
1043:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1044:       };
1045: #endif
1046:       // Specialize to an IConstantSection
1047:       template<typename SendOverlap, typename RecvOverlap, typename Value>
1048:       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) {
1049:         // Why doesn't this work?
1050:         //   This supposed to be a copy, BUT filtered through the sendOverlap
1051:         //   However, an IConstant section does not allow filtration of its
1052:         //   chart. Therefore, you end up with either
1053:         //
1054:         //   a) Too many items in the chart, copied from the remote sendSection
1055:         //   b) A chart mapped to the local numbering, which we do not want
1056:         copyIConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1057:       };
1058:       // Specialize to an BaseSection/ConstantSection pair
1059: #if 0
1060:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
1061:       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) {
1062:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1063:       };
1064: #else
1065:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
1066:       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) {
1067:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1068:       };
1069: #endif
1070:       // Specialize to an BaseSectionV/ConstantSection pair
1071: #if 0
1072:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
1073:       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) {
1074:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1075:       };
1076: #else
1077:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_>
1078:       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) {
1079:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1080:       };
1081: #endif
1082:       // Specialize to an LabelBaseSection/ConstantSection pair
1083: #if 0
1084:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_, typename Label_>
1085:       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) {
1086:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1087:       };
1088: #else
1089:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_, typename Label_>
1090:       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) {
1091:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1092:       };
1093: #endif
1094:       template<typename SendOverlap, typename RecvOverlap, typename Sieve_, typename Label_>
1095:       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) {
1096:         copyConstant(sendOverlap, recvOverlap, sendSection, recvSection);
1097:       };
1098:       // Specialize to a ConstantSection for ArrowSection
1099:       template<typename SendOverlap, typename RecvOverlap, typename Value>
1100:       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) {
1101:         copyConstantArrow(sendOverlap, recvOverlap, sendSection, recvSection);
1102:       };
1103:     };
1104:     class BinaryFusion {
1105:     public:
1106:       template<typename OverlapSection, typename RecvOverlap, typename Section, typename Function>
1107:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section, Function binaryOp) {
1108:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
1109:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

1111:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1112:           // TODO: This must become a more general iterator over colors
1113:           const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
1114:           // Just taking the first value
1115:           const typename Section::point_type&        localPoint    = *p_iter;
1116:           const typename OverlapSection::point_type& remotePoint   = points->begin().color();
1117:           const typename OverlapSection::value_type *overlapValues = overlapSection->restrictPoint(remotePoint);
1118:           const typename Section::value_type        *localValues   = section->restrictPoint(localPoint);
1119:           const int                                  dim           = section->getFiberDimension(localPoint);
1120:           // TODO: optimize allocation
1121:           typename Section::value_type              *values        = new typename Section::value_type[dim];

1123:           for(int d = 0; d < dim; ++d) {
1124:             values[d] = binaryOp(overlapValues[d], localValues[d]);
1125:           }
1126:           section->updatePoint(localPoint, values);
1127:           delete [] values;
1128:         }
1129:       };
1130:     };
1131:     class ReplacementBinaryFusion {
1132:     public:
1133:       template<typename OverlapSection, typename RecvOverlap, typename Section>
1134:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
1135:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
1136:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

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

1145:           section->update(localPoint, overlapSection->restrictPoint(remotePoint));
1146:         }
1147:       };
1148:     };
1149:     class AdditiveBinaryFusion {
1150:     public:
1151:       template<typename OverlapSection, typename RecvOverlap, typename Section>
1152:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
1153:         typedef typename Section::point_type        point_type;
1154:         typedef typename OverlapSection::point_type overlap_point_type;
1155:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
1156:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

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

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

1168:             section->updateAddPoint(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
1169:           }
1170:         }
1171:       };
1172:     };
1173:     class InsertionBinaryFusion {
1174:     public:
1175:       // Insert the overlapSection values into section along recvOverlap
1176:       template<typename OverlapSection, typename RecvOverlap, typename Section>
1177:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<Section>& section) {
1178:         typedef typename Section::point_type        point_type;
1179:         typedef typename OverlapSection::point_type overlap_point_type;
1180:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
1181:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

1183:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1184:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
1185:           const point_type&                              localPoint  = *p_iter;
1186:           const int                                      rank        = *points->begin();
1187:           const point_type&                              remotePoint = points->begin().color();

1189:           if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1190:             if (!section->hasPoint(localPoint)) {
1191:               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;
1192:             }
1193:             section->setFiberDimension(localPoint, overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint)));
1194:           }
1195:         }
1196:         if (rPoints->size()) {section->allocatePoint();}
1197:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1198:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
1199:           const point_type&                              localPoint  = *p_iter;
1200:           const int                                      rank        = *points->begin();
1201:           const point_type&                              remotePoint = points->begin().color();

1203:           if (overlapSection->hasPoint(overlap_point_type(rank, remotePoint))) {
1204:             section->updatePoint(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
1205:           }
1206:         }
1207:       };
1208:       // Specialize to ArrowSection
1209:       template<typename OverlapSection, typename RecvOverlap>
1210:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<UniformSection<MinimalArrow<int,int>,int> >& section) {
1211:         typedef UniformSection<MinimalArrow<int,int>,int> Section;
1212:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
1213:         const typename RecvOverlap::traits::baseSequence::iterator rBegin  = rPoints->begin();
1214:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

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

1221:           for(typename RecvOverlap::traits::baseSequence::iterator q_iter = rBegin; q_iter != rEnd; ++q_iter) {
1222:             const Obj<typename RecvOverlap::coneSequence>& targets      = recvOverlap->cone(*q_iter);
1223:             const typename RecvOverlap::target_type        localTarget  = *q_iter;
1224:             const typename RecvOverlap::target_type        remoteTarget = targets->begin().color();
1225:             const typename Section::point_type             localPoint(localSource, localTarget);
1226:             const typename OverlapSection::point_type      remotePoint(remoteSource, remoteTarget);

1228:             if (overlapSection->hasPoint(remotePoint)) {section->setFiberDimension(localPoint, overlapSection->getFiberDimension(remotePoint));}
1229:           }
1230:         }
1231:         if (rPoints->size()) {section->allocatePoint();}
1232:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rBegin; p_iter != rEnd; ++p_iter) {
1233:           const Obj<typename RecvOverlap::coneSequence>& sources      = recvOverlap->cone(*p_iter);
1234:           const typename RecvOverlap::target_type        localSource  = *p_iter;
1235:           const typename RecvOverlap::target_type        remoteSource = sources->begin().color();

1237:           for(typename RecvOverlap::traits::baseSequence::iterator q_iter = rBegin; q_iter != rEnd; ++q_iter) {
1238:             const Obj<typename RecvOverlap::coneSequence>& targets      = recvOverlap->cone(*q_iter);
1239:             const typename RecvOverlap::target_type        localTarget  = *q_iter;
1240:             const typename RecvOverlap::target_type        remoteTarget = targets->begin().color();
1241:             const typename Section::point_type             localPoint(localSource, localTarget);
1242:             const typename OverlapSection::point_type      remotePoint(remoteSource, remoteTarget);
1243: 
1244:             if (overlapSection->hasPoint(remotePoint)) {section->updatePoint(localPoint, overlapSection->restrictPoint(remotePoint));}
1245:           }
1246:         }
1247:       };
1248:       // Specialize to the Sieve
1249:       template<typename OverlapSection, typename RecvOverlap, typename Renumbering, typename Point>
1250:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, Renumbering& renumbering, const Obj<Sieve<Point,Point,int> >& sieve) {
1251:         typedef typename OverlapSection::point_type overlap_point_type;
1252:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
1253:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

1255:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1256:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
1257:           const Point&                                   localPoint  = *p_iter;
1258:           const int                                      rank        = *points->begin();
1259:           const Point&                                   remotePoint = points->begin().color();
1260:           const int                                      size        = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1261:           const typename OverlapSection::value_type     *values      = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));
1262:           int                                            c           = 0;

1264:           sieve->clearCone(localPoint);
1265:           for(int i = 0; i < size; ++i, ++c) {sieve->addCone(renumbering[values[i]], localPoint, c);}
1266:         }
1267:       };
1268:       // Specialize to the ISieve
1269:       template<typename OverlapSection, typename RecvOverlap, typename Renumbering, typename Point>
1270:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, Renumbering& renumbering, const Obj<IFSieve<Point> >& sieve) {
1271:         typedef typename OverlapSection::point_type overlap_point_type;
1272:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
1273:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();
1274:         int                                                        maxSize = 0;

1276:         for(typename RecvOverlap::traits::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(renumbering[values[i]], 1);}
1286:           for(int i = 0; i < size; ++i) {sieve->addSupportSize(renumbering[values[i].first], 1);}
1287:           maxSize = std::max(maxSize, size);
1288:         }
1289:         sieve->allocate();
1290:         ///typename OverlapSection::value_type *localValues = new typename OverlapSection::value_type[maxSize];
1291:         typename OverlapSection::value_type::first_type  *localValues      = new typename OverlapSection::value_type::first_type[maxSize];
1292:         typename OverlapSection::value_type::second_type *localOrientation = new typename OverlapSection::value_type::second_type[maxSize];

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

1302:           ///for(int i = 0; i < size; ++i) {localValues[i] = renumbering[values[i]];}
1303:           for(int i = 0; i < size; ++i) {
1304:             localValues[i]      = renumbering[values[i].first];
1305:             localOrientation[i] = values[i].second;
1306:           }
1307:           sieve->setCone(localValues, localPoint);
1308:           sieve->setConeOrientation(localOrientation, localPoint);
1309:         }
1310:         delete [] localValues;
1311:         delete [] localOrientation;
1312:       };
1313:       template<typename OverlapSection, typename RecvOverlap, typename Point>
1314:       static void fuse(const Obj<OverlapSection>& overlapSection, const Obj<RecvOverlap>& recvOverlap, const Obj<IFSieve<Point> >& sieve) {
1315:         typedef typename OverlapSection::point_type overlap_point_type;
1316:         const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
1317:         const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();
1318:         int                                                        maxSize = 0;

1320:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1321:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
1322:           const Point&                                   localPoint  = *p_iter;
1323:           const int                                      rank        = *points->begin();
1324:           const Point&                                   remotePoint = points->begin().color();
1325:           const int                                      size        = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1326:           const typename OverlapSection::value_type     *values      = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));

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

1336:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1337:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
1338:           const Point&                                   localPoint  = *p_iter;
1339:           const int                                      rank        = *points->begin();
1340:           const Point&                                   remotePoint = points->begin().color();
1341:           const int                                      size        = overlapSection->getFiberDimension(overlap_point_type(rank, remotePoint));
1342:           const typename OverlapSection::value_type     *values      = overlapSection->restrictPoint(overlap_point_type(rank, remotePoint));

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

1361:         for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1362:           const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
1363:           const typename Section::point_type&            localPoint  = *p_iter;
1364:           const int                                      rank        = *points->begin();
1365:           const typename OverlapSection::point_type&     remotePoint = points->begin().color();

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

1376:           section->update(localPoint, overlapSection->restrictPoint(overlap_point_type(rank, remotePoint)));
1377:         }
1378:       };
1379:     };
1380:   }
1381: }

1383: #endif