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