Actual source code: CoSifter.hh

  1: #ifndef included_ALE_CoSifter_hh
  2: #define included_ALE_CoSifter_hh

  4: #ifndef  included_ALE_Sieve_hh
  5: #include <Sieve.hh>
  6: #endif
  7: #ifndef  included_ALE_ParDelta_hh
  8: #include <ParDelta.hh>
  9: #endif
 10: #include <list>
 11: #include <stack>
 12: #include <queue>

 14: // Dmitry's explanation:
 15: //
 16: // Okay, check out what I have put there.
 17: // It's a rather high-level interface, but I think it sketches out the implementation idea.  I have also become a master of switching from 'public' to 'private' and back.

 19: // The idea is to put more power into Sifters (bipartite graphs with color).  They are like Sieves but with two point types (source and target) and no recursive operations (nCone, closure, etc).
 20: // I claim they should be parallel, so cone/support completions should be computable for them.  The footprint is incorporated into the color of the new Sifter, which is returned as a completion.
 21: // It would be very natural to have Sieve<Point_, Color_> to extend Sifter<Point_, Point_, Color_> with the recursive operations.

 23: // The reason for putting the completion functionality into Sifters is that patches and indices under and over a topology Sieve are Sifters and have to be completed:
 24: // the new overlap_patches has to encode patch pairs along with the rank of the second patch (first is always local); likewise, overlap_indices must encode a pair of intervals with a rank
 25: // -- the attached to the same Sieve point by two different processes -- one local and one (possibly) remote.  At any rate, the support completion of 'patches' contains all the information
 26: // needed for 'overlap_patches' -- remember that struct with a triple {point, patch, number} you had on the board?   Likewise for 'overlap_indices' built out of the cone completion of 'indices'.

 28: // Once the 'overlap_XXX' are computed, we can allocate the storage for the Delta data and post sends receives.
 29: // We should be able to reuse the completion subroutine from the old Sieve.
 30: // So you are right that perhaps Sieve completion gets us to the CoSieve completion, except I think it's the Sifter completion this time.
 31: // I can do the completion when I come back if you get the serial Sifter/Sieve stuff going.
 32: //
 33: namespace ALE {
 34:     template <typename Sieve_, typename Patch_, typename Index_, typename Value_>
 35:     class CoSifter {
 36:       template <class _Type>
 37:       struct trueFunc : public std::binary_function<_Type, _Type, bool>
 38:       {
 39:         bool operator()(const _Type& __x, const _Type& __y) const { return true; }
 40:       };
 41:     public:
 42:       // Basic types
 43:       typedef Sieve_ sieve_type;
 44:       typedef typename sieve_type::point_type point_type;
 45:       typedef std::vector<point_type> PointArray;
 46:       typedef Patch_ patch_type;
 47:       typedef Index_ index_type;
 48:       typedef std::vector<index_type> IndexArray;
 49:       typedef Value_ value_type;
 50:       typedef Sifter<point_type,patch_type,index_type, ::boost::multi_index::composite_key_compare<std::less<point_type>, trueFunc<index_type>, std::less<patch_type> >, SifterDef::RecContainer<point_type, SifterDef::Rec<point_type> >, SifterDef::RecContainer<patch_type, SifterDef::Rec<patch_type> > > order_type;

 52:       typedef RightSequenceDuplicator<ConeArraySequence<typename sieve_type::traits::arrow_type> > fuser;
 53:       typedef ParConeDelta<sieve_type, fuser,
 54:                            typename sieve_type::template rebind<typename fuser::fusion_source_type,
 55:                                                                 typename fuser::fusion_target_type,
 56:                                                                 typename fuser::fusion_color_type,
 57:                                                                 typename sieve_type::traits::cap_container_type::template rebind<typename fuser::fusion_source_type,
 58:                                                                                                                                  typename sieve_type::traits::sourceRec_type::template rebind<typename fuser::fusion_source_type,
 59:                                                                                                                                                                                               typename sieve_type::marker_type>::type>::type,
 60:                                                                 typename sieve_type::traits::base_container_type::template rebind<typename fuser::fusion_target_type,
 61:                                                                                                                                   typename sieve_type::traits::targetRec_type::template rebind<typename fuser::fusion_target_type,
 62:                                                                                                                                                                                                typename sieve_type::marker_type>::type>::type>::type> coneDelta_type;
 63:       typedef ParSupportDelta<sieve_type, fuser,
 64:                               typename sieve_type::template rebind<typename fuser::fusion_source_type,
 65:                                                                    typename fuser::fusion_target_type,
 66:                                                                    typename fuser::fusion_color_type,
 67:                                                                    typename sieve_type::traits::cap_container_type::template rebind<typename fuser::fusion_source_type, typename sieve_type::traits::sourceRec_type::template rebind<typename fuser::fusion_source_type, typename sieve_type::marker_type>::type>::type,
 68:                                                                    typename sieve_type::traits::base_container_type::template rebind<typename fuser::fusion_target_type, typename sieve_type::traits::targetRec_type::template rebind<typename fuser::fusion_target_type, typename sieve_type::marker_type>::type>::type
 69:       >::type> supportDelta_type;
 70:       typedef RightSequenceDuplicator<ConeArraySequence<typename order_type::traits::arrow_type> > orderFuser;
 71:       typedef ParSupportDelta<order_type, orderFuser> supportOrderDelta_type;
 72:       typedef CoSifter<sieve_type, patch_type, index_type, int> bundle_type;
 73:     private:
 74:       MPI_Comm        _comm;
 75:       int             _commRank;
 76:       int             _commSize;
 77:       Obj<sieve_type> _topology;
 78:       // We need an ordering, which should be patch<--order--point
 79:       Obj<order_type> _order;
 80:       // We need a reordering, which should be patch<--new order--old order
 81:       std::map<std::string,Obj<order_type> > _reorders;
 82:       // We can add fields to an ordering using <patch,field><--order--point
 83:       // We need sequences that can return the color, or do it automatically
 84:       // We allocate based upon a certain
 85:       std::map<patch_type,int>          _storageSize;
 86:       std::map<patch_type,value_type *> _storage;
 87:       int *offsets;
 88:       int  ghostVars;
 89:       Obj<bundle_type> localOrder;
 90:       Obj<bundle_type> globalOrder;
 91:     public:
 92:       int              debug;
 93:       // OLD CRAP:
 94:       // Breakdown of the base Sieve into patches
 95:       //   the colors (int) order the points (point_type) over a patch (patch_type).
 96:       // A patch is a member of the sheaf over the sieve which indicates a maximal
 97:       // domain of definition for a function on the sieve. Furthermore, we use the
 98:       // patch coloring to order the patch values, the analog of a coordinate system
 99:       // on the patch. We use a Sifter here, but the object can properly be thought
100:       // of as a CoSieve over the topology sieve.
101:     public:
102:       CoSifter(MPI_Comm comm = PETSC_COMM_SELF, int debug = 0) : _comm(comm), debug(debug) {
103:         this->_order = new order_type(this->_comm, debug);
104:         MPI_Comm_rank(this->_comm, &this->_commRank);
105:         MPI_Comm_size(this->_comm, &this->_commSize);
106:         this->offsets = NULL;
107:         this->ghostVars = 0;
108:       };
109:       ~CoSifter() {
110:         if (this->offsets) {
111:           delete [] this->offsets;
112:           this->offsets = NULL;
113:         }
114:       };

116:       MPI_Comm        comm() const {return this->_comm;};
117:       int             commRank() const {return this->_commRank;};
118:       int             commSize() const {return this->_commSize;};
119:       void            setTopology(const Obj<sieve_type>& topology) {this->_topology = topology;};
120:       Obj<sieve_type> getTopology() const {return this->_topology;};
121:       // -- Patch manipulation --
122:       // Creates a patch whose order is taken from the input point sequence
123:       template<typename pointSequence> void setPatch(const Obj<pointSequence>& points, const patch_type& patch) {
124:         int c = 1;

126:         for(typename pointSequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
127:           this->_order->addArrow(*p_iter, patch, index_type(c++, 0));
128:         }
129:         if (points->begin() == points->end()) {
130:           this->_order->addBasePoint(patch);
131:         }
132:       };
133:       void extendPatch(const point_type& point, const patch_type& patch, const index_type& color) {
134:         this->_order->addArrow(point, patch, color);
135:       };
136:       void extendPatch(const point_type& point, const patch_type& patch) {
137:         Obj<typename order_type::traits::coneSequence> cone = this->_order->cone(patch);
138:         int c = 1;

140:         for(typename order_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
141:           c = c_iter.color().prefix;
142:         }
143:         this->extendPatch(point, patch, point_type(c+1, 0));
144:       };
145:       Obj<bundle_type> getGlobalOrder() const {
146:         return this->globalOrder;
147:       };
148:       Obj<bundle_type> getLocalOrder() const {
149:         return this->localOrder;
150:       };
151:     public:
152:       Obj<order_type> __getOrder() {
153:         return this->_order;
154:       };
155:       Obj<order_type> __getOrder(const std::string& orderName) {
156:         if (this->_reorders.find(orderName) == this->_reorders.end()) {
157:           if (this->debug) {std::cout << "Creating new order: " << orderName << std::endl;}
158:           this->_reorders[orderName] = new order_type(this->_comm, this->debug);
159:         }
160:         return this->_reorders[orderName];
161:       };
162:     public:
163:       // Creates a patch for a named reordering whose order is taken from the input point sequence
164:       void setPatch(const std::string& orderName, const point_type& point, const patch_type& patch) {
165:         Obj<order_type> reorder = this->__getOrder(orderName);

167:         reorder->addArrow(point, patch, point_type());
168:       }
169:       template<typename pointSequence> void setPatch(const std::string& orderName, const Obj<pointSequence>& points, const patch_type& patch) {
170:         Obj<order_type> reorder = this->__getOrder(orderName);
171:         int c = 1;

173:         for(typename pointSequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
174:           reorder->addArrow(*p_iter, patch, point_type(c++, 0));
175:         }
176:       };
177:       // Returns the points in the patch in order
178:       Obj<typename order_type::coneSequence> getPatch(const patch_type& patch) const {
179:         return this->_order->cone(patch);
180:       };
181:       Obj<typename order_type::baseSequence> getPatches() {
182:         return this->_order->base();
183:       };
184:     private:
185:       void __checkOrderName(const std::string& orderName) {
186:         if (this->_reorders.find(orderName) != this->_reorders.end()) return;
187:         //FIX: String memory management
188:         std::string msg("Invalid order name: ");

190:         msg += orderName;
191:         throw ALE::Exception(msg.c_str());
192:       };
193:     public:
194:       // Returns the points in the reorder patch in order
195:       Obj<typename order_type::coneSequence> getPatch(const std::string& orderName, const patch_type& patch) {
196:         this->__checkOrderName(orderName);
197:         return this->_reorders[orderName]->cone(patch);
198:       };
199:       Obj<typename order_type::baseSequence> getPatches(const std::string& orderName) {
200:         return this->_reorders[orderName]->base();
201:       };
202:       // -- Index manipulation --
203:     private:
204:       struct changeOffset {
205:         changeOffset(int newOffset) : newOffset(newOffset) {};

207:         void operator()(typename order_type::Arrow_& p) const {
208:           p.color.prefix = newOffset;
209:         }
210:       private:
211:         int newOffset;
212:       };
213:       struct incrementOffset {
214:         incrementOffset(int newOffset) : newOffset(newOffset) {};

216:         void operator()(typename order_type::Arrow_& p) const {
217:           p.color.prefix = p.color.prefix + newOffset;
218:         }
219:       private:
220:         int newOffset;
221:       };
222:       struct changeDim {
223:         changeDim(int newDim) : newDim(newDim) {};

225:         void operator()(typename order_type::Arrow_& p) const {
226:           p.color.index = newDim;
227:         }
228:       private:
229:         int newDim;
230:       };
231:       struct changeIndex {
232:         changeIndex(int newOffset, int newDim) : newOffset(newOffset), newDim(newDim) {};

234:         void operator()(typename order_type::Arrow_& p) const {
235:           p.color.prefix = newOffset;
236:           p.color.index  = newDim;
237:         }
238:       private:
239:         int newOffset;
240:         int newDim;
241:       };
242:     public:
243:       int getFiberDimension(const patch_type& patch, const point_type& p) const {
244:         return this->_order->getColor(p, patch, false).index;
245:       };
246:       int getFiberDimension(const std::string& orderName, const patch_type& patch, const point_type& p) const {
247:         this->__checkOrderName(orderName);
248:         return this->_reorders[orderName]->getColor(p, patch, false).index;
249:       };
250:       void setFiberDimension(const patch_type& patch, const point_type& p, int dim) {
251:         this->_order->modifyColor(p, patch, changeDim(-dim));
252:       };
253:       void setFiberDimensionByDepth(const patch_type& patch, int depth, int dim) {
254:         Obj<typename sieve_type::traits::depthSequence> points = this->_topology->depthStratum(depth);

256:         for(typename sieve_type::traits::depthSequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
257:           this->setFiberDimension(patch, *p_iter, dim);
258:         }
259:       };
260:       void setFiberDimensionByHeight(const patch_type& patch, int height, int dim) {
261:         Obj<typename sieve_type::traits::heightSequence> points = this->_topology->heightStratum(height);

263:         for(typename sieve_type::traits::heightSequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
264:           this->setFiberDimension(patch, *p_iter, dim);
265:         }
266:       };
267:       void setFiberDimension(const std::string& orderName, const patch_type& patch, const point_type& p, int dim) {
268:         this->__checkOrderName(orderName);
269:         this->_reorders[orderName]->modifyColor(p, patch, changeDim(-dim));
270:       };
271:       void setFiberDimensionByDepth(const std::string& orderName, const patch_type& patch, int depth, int dim) {
272:         Obj<typename sieve_type::depthSequence> points = this->_topology->depthStratum(depth);

274:         for(typename sieve_type::depthSequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
275:           this->setFiberDimension(orderName, patch, *p_iter, dim);
276:         }
277:       };
278:       int getFiberOffset(const patch_type& patch, const point_type& p) const {
279:         return this->_order->getColor(p, patch, false).prefix;
280:       };
281:       void setFiberOffset(const patch_type& patch, const point_type& p, int offset) {
282:         this->_order->modifyColor(p, patch, changeOffset(offset));
283:       };
284:       void addFiberOffset(const patch_type& patch, const point_type& p, int offset) {
285:         this->_order->modifyColor(p, patch, incrementOffset(offset));
286:       };
287:     private:
288:       struct trueTester {
289:       public:
290:         bool operator()(const point_type& p) const {
291:           return true;
292:         };
293:       };
294:       template<typename OrderTest>
295:       void __orderCell(const Obj<order_type>& order, const patch_type& patch, const point_type& cell, int& offset, const OrderTest& tester) {
296:         // Set the prefix to the current offset (this won't kill the topology iterator)
297:         Obj<typename sieve_type::coneSequence> cone = this->_topology->cone(cell);
298:         typename sieve_type::coneSequence::iterator end = cone->end();

300:         for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
301:           this->__orderCell(order, patch, *p_iter, offset, tester);
302:         }

304:         int dim = order->getColor(cell, patch, false).index;

306:         if ((dim < 0) && tester(cell)) {
307:           order->modifyColor(cell, patch, changeIndex(offset, -dim));
308:           if (debug) {std::cout << "Order point " << cell << " of size " << -dim << " and offset " << offset << " color " << order->getColor(cell, patch) << std::endl;}
309:           offset -= dim;
310:         }
311:       };
312:       // This constructs an order on the patch by fusing the Ord CoSieve (embodied by the prefix number)
313:       // and the Count CoSieve (embodied by the index), turning the prefix into an offset.
314:       template<typename OrderTest>
315:       void __orderPatch(const Obj<order_type>& order, const patch_type& patch, bool allocate, const OrderTest& tester, const PointArray& points) {
316:         int offset = 0;

318:         // Loop over patch members
319:         for(typename PointArray::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
320:           // Traverse the closure of the member in the topology
321:           if (debug) {std::cout << "Ordering patch point " << *p_iter << std::endl;}
322:           this->__orderCell(order, patch, *p_iter, offset, tester);
323:         }
324:         if (allocate) {
325:           this->setPatchSize(patch, offset);
326:         }
327:       };
328:       // Filter out newly added points and then number
329:       template<typename OrderTest>
330:       void __orderPatch(const Obj<order_type>& order, const patch_type& patch, bool allocate, const OrderTest& tester) {
331:         PointArray points;

333:         Obj<typename order_type::coneSequence> cone = order->cone(patch);
334:         int rank = 1;
335:         for(typename order_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
336:           if (p_iter.color().prefix == rank) {
337:             points.push_back(*p_iter);
338:             rank++;
339: //           } else if (debug) {
340: //             std::cout << "Rejected patch point " << *p_iter << " with color " << p_iter.color() << std::endl;
341:           }
342:         }
343:         this->__orderPatch(order, patch, allocate, tester, points);
344:       };
345:       // Number all points
346:       template<typename OrderTest>
347:       void __reorderPatch(const Obj<order_type>& order, const patch_type& patch, bool allocate, const OrderTest& tester) {
348:         PointArray points;

350:         Obj<typename order_type::coneSequence> cone = order->cone(patch);
351:         for(typename order_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
352:           points.push_back(*p_iter);
353:         }
354:         int rank = 1;
355:         for(typename PointArray::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
356:           order->modifyColor(*p_iter, patch, changeIndex(rank++, -(order->getColor(*p_iter, patch, false).index)));
357:         }
358:         this->__orderPatch(order, patch, allocate, tester, points);
359:       };
360:     public:
361:       void setPatchSize(const patch_type& patch, int size = -1) {
362:         if (size < 0) {
363:           Obj<typename order_type::coneSequence> cone = getPatch(patch);

365:           for(typename order_type::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
366:             const index_type& idx = this->getIndex(patch, *c_iter);

368:             if (size < idx.prefix + idx.index) size = idx.prefix + idx.index;
369:           }
370:         }
371:         this->_storageSize[patch] = size;
372:       };
373:       void allocatePatches() {
374:         Obj<typename order_type::baseSequence> patches = this->getPatches();
375:         int size = 0;

377:         for(typename order_type::baseSequence::iterator p_iter = patches->begin(); p_iter != patches->end(); ++p_iter) {
378:           size += this->_storageSize[*p_iter];
379:         }
380:         if (this->_storage.find(*patches->begin()) != this->_storage.end()) {
381:           delete [] this->_storage[*patches->begin()];
382:         }
383:         value_type *data = new value_type[size];
384:         memset(data, 0, size*sizeof(value_type));
385:         for(typename order_type::baseSequence::iterator p_iter = patches->begin(); p_iter != patches->end(); ++p_iter) {
386:           this->_storage[*p_iter] = data;
387:           data += this->_storageSize[*p_iter];
388:         }
389:       };
390:       void orderPatch(const patch_type& patch) {
391:         this->__orderPatch(this->_order, patch, true, trueTester());
392:       }
395:       template<typename OrderTest>
396:       void orderPatches(const OrderTest& tester) {
397:         ALE_LOG_EVENT_BEGIN;
398:         Obj<typename order_type::baseSequence> base = this->_order->base();

400:         for(typename order_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
401:           this->__orderPatch(this->_order, *b_iter, true, tester);
402:         }
403:         this->allocatePatches();
404:         ALE_LOG_EVENT_END;
405:       };
406:       void orderPatches() {
407:         this->orderPatches(trueTester());
408:       };
409:       void orderPatches(const std::string& orderName) {
410:         this->__checkOrderName(orderName);
411:         ALE_LOG_EVENT_BEGIN;
412:         Obj<typename order_type::baseSequence> base = this->_reorders[orderName]->base();

414:         for(typename order_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
415:           this->__orderPatch(this->_reorders[orderName], *b_iter, false, trueTester());
416:         }
417: //         std::cout << orderName << " ordering:" << std::endl;
418: //         for(typename order_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
419: //           Obj<typename order_type::coneSequence> cone = this->getPatch(orderName, *b_iter);

421: //           std::cout << "  patch " << *b_iter << std::endl;
422: //           for(typename order_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
423: //             std::cout << "    " << *p_iter << std::endl;
424: //           }
425: //         };
426:         ALE_LOG_EVENT_END;
427:       };
430:       template<typename OrderTest>
431:       void reorderPatches(const OrderTest& tester) {
432:         ALE_LOG_EVENT_BEGIN;
433:         Obj<typename order_type::baseSequence> base = this->_order->base();

435:         for(typename order_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
436:           this->__reorderPatch(this->_order, *b_iter, true, tester);
437:         }
438:         this->allocatePatches();
439:         ALE_LOG_EVENT_END;
440:       };
441:       void reorderPatches() {
442:         this->reorderPatches(trueTester());
443:       };
444:       //Obj<IndexArray> getIndices(const patch_type& patch);
445:       //Obj<IndexArray> getIndices(const patch_type& patch, const point_type& p);
446:       const index_type& getIndex(const patch_type& patch, const point_type& p) {
447:         return this->_order->getColor(p, patch, false);
448:       };
449:       Obj<IndexArray> getIndices(const std::string& orderName, const patch_type& patch) {
450:         Obj<typename order_type::coneSequence> cone = getPatch(orderName, patch);
451:         Obj<IndexArray>                        array = new IndexArray();
452:         patch_type                             oldPatch;

454:         // We have no way to map the the old patch yet
455:         // It would be better to map this through in a sequence to the original indices (like fusion)
456:         for(typename order_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
457:           array->push_back(this->getIndex(oldPatch, *p_iter));
458:         }
459:         return array;
460:       }
461:       // -- Value manipulation --
462:     private:
463:       void __checkPatch(const patch_type& patch) const {
464:         if (this->_storage.find(patch) != this->_storage.end()) return;
465:         ostringstream msg;

467:         msg << "Invalid patch: " << patch;
468:         throw ALE::Exception(msg.str().c_str());
469:       };
470:     public:
471:       int getSize() {
472:         Obj<typename order_type::baseSequence> patches = this->getPatches();
473:         int totalSize = 0;

475:         for(typename order_type::baseSequence::iterator p_iter = patches->begin(); p_iter != patches->end(); ++p_iter) {
476:           totalSize += this->_storageSize[*p_iter];
477:         }
478:         return totalSize;
479:       };
480:       int getSize(const patch_type& patch) {
481:         this->__checkPatch(patch);
482:         return this->_storageSize[patch];
483:       };
484:       const int *getGlobalOffsets() {return this->offsets;};
485:       const value_type *restrict(const patch_type& patch, bool fail = true) {
486:         if (fail) {
487:           this->__checkPatch(patch);
488:         } else {
489:           if (this->_storage.find(patch) == this->_storage.end()) {
490:             return NULL;
491:           }
492:         }
493:         return this->_storage[patch];
494:       };
495:       const value_type *restrict(const patch_type& patch, const point_type& p) {
496:         this->__checkPatch(patch);
497:         return &this->_storage[patch][this->_order->getColor(p, patch, false).prefix];
498:       };
499:       // Can this be improved?
500:       const value_type *restrict(const std::string& orderName, const patch_type& patch) {
501:         Obj<typename order_type::coneSequence> cone = getPatch(orderName, patch);
502:         static value_type                     *values = NULL;
503:         static int                             size = 0;
504:         int                                    newSize = 0;
505:         int                                    newI = 0;
506:         patch_type                             oldPatch;

508:         for(typename order_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
509:           newSize += this->getIndex(oldPatch, *p_iter).index;
510:         }
511:         if (newSize != size) {
512:           if (!values) delete [] values;
513:           size = newSize;
514:           values = new value_type[size];
515:         }
516:         for(typename order_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
517:           const index_type& ind = this->getIndex(oldPatch, *p_iter);

519:           for(int i = ind.prefix; i < ind.prefix+ind.index; ++i) {
520:             values[newI++] = this->_storage[oldPatch][i];
521:           }
522:         }
523:         return values;
524:       };
525:       const value_type *restrict(const std::string& orderName, const patch_type& patch, const point_type& p);
526:       void              update(const patch_type& patch, const value_type values[]) {
527:         this->__checkPatch(patch);
528:         value_type *storage = this->_storage[patch];
529:         const int   size = this->_storageSize[patch];

531:         for(int i = 0; i < size; ++i) {
532:           storage[i] = values[i];
533:         }
534:       };
535:       void              update(const patch_type& patch, const point_type& p, const value_type values[]) {
536:         const index_type& idx = this->getIndex(patch, p);
537:         int offset = idx.prefix;
538:         value_type *storage = &(this->_storage[patch][offset]);

540:         for(int i = 0; i < idx.index; ++i) {
541:           if (debug) {std::cout << "Set a[" << offset+i << "] = " << values[i] << " on patch " << patch << std::endl;}
542:           storage[i] = values[i];
543:         }
544:       };
545:       // Can this be improved?
546:       void              update(const std::string& orderName, const patch_type& patch, const value_type values[]) {
547:         Obj<typename order_type::coneSequence> cone = getPatch(orderName, patch);
548:         int                                    newI = 0;
549:         patch_type                             oldPatch;

551:         for(typename order_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
552:           const index_type& ind = this->getIndex(oldPatch, *p_iter);

554:           for(int i = ind.prefix; i < ind.prefix+ind.index; ++i) {
555:             this->_storage[oldPatch][i] = values[newI++];
556:           }
557:         }
558:       };
559:       void              update(const std::string& orderName, const patch_type& patch, const point_type& p, const value_type values[]);
560:       void              updateAdd(const patch_type& patch, const value_type values[]);
561:       void              updateAdd(const patch_type& patch, const point_type& p, const value_type values[]) {
562:         const index_type& idx = this->getIndex(patch, p);
563:         int offset = idx.prefix;

565:         for(int i = 0; i < idx.index; ++i) {
566:           if (debug) {std::cout << "Set a[" << offset+i << "] = " << values[i] << " on patch " << patch << std::endl;}
567:           this->_storage[patch][offset+i] += values[i];
568:         }
569:       };
570:       void              updateAdd(const std::string& orderName, const patch_type& patch, const value_type values[]) {
571:         Obj<typename order_type::coneSequence> cone = getPatch(orderName, patch);
572:         int                                    newI = 0;
573:         patch_type                             oldPatch;

575:         for(typename order_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
576:           const index_type& ind = this->getIndex(oldPatch, *p_iter);

578:           for(int i = ind.prefix; i < ind.prefix+ind.index; ++i) {
579:             this->_storage[oldPatch][i] += values[newI++];
580:           }
581:         }
582:       };
583:       void              updateAdd(const std::string& orderName, const patch_type& patch, const point_type& p, const value_type values[]);

585:       void view(const char* label) const {
586:         ostringstream txt;

588:         if(label != NULL) {
589:           if(this->commRank() == 0) {
590:             txt << "viewing CoSifter :'" << label << "'" << std::endl;
591:           }
592:         } else {
593:           if(this->commRank() == 0) {
594:             txt << "viewing a CoSifter" << std::endl;
595:           }
596:         }
597:         for(typename std::map<patch_type,value_type *>::const_iterator s_iter = this->_storage.begin(); s_iter != this->_storage.end(); ++s_iter) {
598:           patch_type patch = s_iter->first;

600:           txt << "[" << this->commRank() << "]: Patch " << patch << std::endl;
601:           Obj<typename order_type::coneSequence> cone = this->getPatch(s_iter->first);
602:           const value_type *array = ((std::map<patch_type,value_type *>) _storage)[s_iter->first];

604:           for(typename order_type::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
605:             index_type color = this->_order->getColor(*c_iter, s_iter->first, false);

607:             if (color.index != 0) {
608:               txt << "[" << this->commRank() << "]:   " << *c_iter << " dim " << color.index << " offset " << color.prefix << "  ";
609:               for(int i = 0; i < color.index; i++) {
610:                 txt << " " << array[color.prefix+i];
611:               }
612:               txt << std::endl;
613:             }
614:           }
615:         }
616:         PetscSynchronizedPrintf(this->comm(), txt.str().c_str());
617:         PetscSynchronizedFlush(this->comm());
618:       };
619:     protected:
620:       struct supportLocalizer {
621:         typedef typename supportDelta_type::overlap_type                      overlap_type;
622:         typedef typename supportDelta_type::overlap_type::traits::capSequence sequence_type;
623:         Obj<overlap_type>  overlap;
624:         Obj<sequence_type> points;
625:         int                rank;
626:         public:
627:         supportLocalizer(Obj<overlap_type> overlap, const int rank) : overlap(overlap), rank(rank) {points = overlap->cap();};
628:         bool isLocal(const Obj<overlap_type> overlap, const Obj<sequence_type> points, const typename sieve_type::point_type& p) const {
629:           if (points->contains(p)) {
630:             Obj<typename overlap_type::traits::supportSequence> neighbors = overlap->support(p);

632:             for(typename overlap_type::traits::supportSequence::iterator s_iter = neighbors->begin(); s_iter != neighbors->end(); ++s_iter) {
633:               if (s_iter.target() < rank)
634:                 return false;
635:             }
636:           }
637:           return true;
638:         };

640:         bool operator()(const typename sieve_type::point_type& p) const {
641:           //std::cout << "Checking for local point " << p << std::endl;
642:           return this->isLocal(this->overlap, this->points, p);
643:         }
644:       };
645:       //FIX: Should just flip internals I think
646:       struct coneLocalizer {
647:         typedef typename coneDelta_type::overlap_type                       overlap_type;
648:         typedef typename coneDelta_type::overlap_type::traits::baseSequence sequence_type;
649:         Obj<overlap_type>  overlap;
650:         Obj<sequence_type> points;
651:         int                rank;
652:         public:
653:         coneLocalizer(Obj<overlap_type> overlap, const int rank) : overlap(overlap), rank(rank) {points = overlap->base();};
654:         bool isLocal(const Obj<overlap_type> overlap, const Obj<sequence_type> points, const typename sieve_type::point_type& p) const {
655:           if (points->contains(p)) {
656:             Obj<typename overlap_type::traits::coneSequence> neighbors = overlap->cone(p);

658:             for(typename overlap_type::traits::coneSequence::iterator c_iter = neighbors->begin(); c_iter != neighbors->end(); ++c_iter) {
659:               if (c_iter.source() < rank)
660:                 return false;
661:             }
662:           }
663:           return true;
664:         };

666:         bool operator()(const typename sieve_type::point_type& p) const {
667:           //std::cout << "Checking for local point " << p << std::endl;
668:           return this->isLocal(this->overlap, this->points, p);
669:         }
670:       };
671:     public:
672:       void createGlobalOrder() {
673:         Obj<typename sieve_type::traits::depthSequence>  vertices = this->_topology->depthStratum(0);
674:         Obj<typename sieve_type::traits::heightSequence> cells    = this->_topology->heightStratum(0);
675:         int useLocBaseOverlap = 0, useLocCapOverlap = 0;
676:         int useBaseOverlap    = 0, useCapOverlap    = 0;
677:         typename bundle_type::patch_type patch;

679:         for(typename sieve_type::traits::depthSequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
680:           Obj<typename order_type::baseSequence> patches = this->getPatches();

682:           for(typename order_type::baseSequence::iterator p_iter = patches->begin(); p_iter != patches->end(); ++p_iter) {
683:             if (this->getFiberDimension(*p_iter, *v_iter) > 0) {
684:               useLocCapOverlap = 1;
685:               break;
686:             }
687:           }
688:           if (useLocCapOverlap) break;
689:         }
690:         for(typename sieve_type::traits::heightSequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
691:           Obj<typename order_type::baseSequence> patches = this->getPatches();

693:           for(typename order_type::baseSequence::iterator p_iter = patches->begin(); p_iter != patches->end(); ++p_iter) {
694:             if (this->getFiberDimension(*p_iter, *c_iter) > 0) {
695:               useLocBaseOverlap = 1;
696:               break;
697:             }
698:           }
699:           if (useLocBaseOverlap) break;
700:         }
701:         // Get overlap
702:         Obj<typename supportDelta_type::overlap_type> capOverlap;
703:         Obj<typename coneDelta_type::overlap_type> baseOverlap;
704:         int rank = this->commRank();

706:         MPI_Allreduce(&useLocCapOverlap, &useCapOverlap, 1, MPI_INT, MPI_LOR, this->comm());
707:         if (useCapOverlap) {
708:           if (this->debug) {
709:             std::cout << "Doing cap overlap" << std::endl;
710:           }
711:           //capOverlap = supportDelta_type::overlap(this->_topology);
712:         }
713:         MPI_Allreduce(&useLocBaseOverlap, &useBaseOverlap, 1, MPI_INT, MPI_LOR, this->comm());
714:         if (useBaseOverlap) {
715:           if (this->debug) {
716:             std::cout << "Doing base overlap" << std::endl;
717:           }
718:           //baseOverlap = coneDelta_type::overlap(this->_topology);
719:         }
720:         if (useCapOverlap && useBaseOverlap) {
721:           throw ALE::Exception("Cannot have both kinds of overlap");
722:         }
723:         if (!useCapOverlap && !useBaseOverlap) {
724:           throw ALE::Exception("Cannot have no overlap");
725:         }
726:         // Give a local offset to each local element, continue sequential offsets for ghosts
727:         // Local order is a CoSifter<sieve_type, patch_type, point_type, int>
728:         //   which means localOrder->_order is a Sifter<point_type,patch_type,point_type>
729:         // SupportDelta::overlap_type is an ASifter<ALE::Point, int, ALE::pair<int,int>, uniColor>
730:         this->localOrder  = bundle_type(this->_comm, this->debug);
731:         this->globalOrder = bundle_type(this->_comm, this->debug);
732:         Obj<typename order_type::baseSequence> base = this->_order->base();

734:         this->localOrder->setTopology(this->_topology);
735:         this->globalOrder->setTopology(this->_topology);
736:         for(typename order_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
737:           Obj<typename order_type::coneSequence> cone = this->getPatch(*b_iter);

739:           this->localOrder->setPatch(cone, *b_iter);
740:           this->globalOrder->setPatch(cone, *b_iter);
741:           for(typename order_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
742:             this->localOrder->setFiberDimension(*b_iter, *p_iter, this->getFiberDimension(*b_iter, *p_iter));
743:             this->globalOrder->setFiberDimension(*b_iter, *p_iter, this->getFiberDimension(*b_iter, *p_iter));
744:           }
745:         }
746: //         if (useCapOverlap) {
747: //           this->localOrder->orderPatches(supportLocalizer(capOverlap, this->commRank()));
748: //           this->globalOrder->orderPatches(supportLocalizer(capOverlap, this->commRank()));
749: //         } else {
750: //           this->localOrder->orderPatches(coneLocalizer(baseOverlap, this->commRank()));
751: //           this->globalOrder->orderPatches(coneLocalizer(baseOverlap, this->commRank()));
752: //         }
753:         if (this->debug) {
754:           this->localOrder->view("Local order");
755:           this->globalOrder->view("Global order");
756:         }
757:         int ierr, localVars = 0;

759:         if (this->offsets) {
760:           delete [] this->offsets;
761:         }
762:         this->offsets = new int[this->_commSize+1];
763:         for(typename order_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
764:           localVars += this->globalOrder->getSize(*b_iter);
765:         }
766:         MPI_Allgather(&localVars, 1, MPI_INT, &this->offsets[1], 1, MPI_INT, this->comm());CHKERROR(ierr, "Error in MPI_Allgather");
767:         this->offsets[0] = 0;
768:         this->ghostVars  = 0;
769:         for(int p = 1; p <= this->commSize(); p++) {
770:           this->offsets[p] += this->offsets[p-1];
771:         }
772:         // Create global numbering
773:         for(typename order_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
774:           Obj<typename order_type::coneSequence> cone = getPatch(*b_iter);

776:           for(typename order_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
777:             this->globalOrder->addFiberOffset(*b_iter, *p_iter, this->offsets[rank]);
778:             this->localOrder->setFiberOffset(*b_iter, *p_iter, this->offsets[rank] + this->ghostVars++);
779:           }
780:         }
781:         if (this->debug) {
782:           this->localOrder->view("Local order with offset");
783:           this->globalOrder->view("Global order with offset");
784:         }
785:         // Complete order to get ghost offsets
786:         //this->globalOrder->completeOrder();
787:         if (this->debug) {
788:           this->globalOrder->view("Global order after completion");
789:         }
790:       };
791:     };
792: } // namespace ALE

794: #endif