Actual source code: Numbering.hh

  1: #ifndef included_ALE_Numbering_hh
  2: #define included_ALE_Numbering_hh

  4: #ifndef  included_ALE_SectionCompletion_hh
  5: #include <SectionCompletion.hh>
  6: #endif


  9: namespace ALE {
 10:     // We have a dichotomy between \emph{types}, describing the structure of objects,
 11:     //   and \emph{concepts}, describing the role these objects play in the algorithm.
 12:     //   Below we identify concepts with potential implementing types.
 13:     //
 14:     //   Concept           Type
 15:     //   -------           ----
 16:     //   Overlap           Sifter
 17:     //   Atlas             ConstantSection, UniformSection
 18:     //   Numbering         UniformSection
 19:     //   GlobalOrder       UniformSection
 20:     //
 21:     // We will use factory types to create objects which satisfy a given concept.
 22:   template<typename Point_, typename Value_ = int, typename Alloc_ = malloc_allocator<Point_> >
 23:   class Numbering : public UniformSection<Point_, Value_, 1, Alloc_> {
 24:   public:
 25:     typedef UniformSection<Point_, Value_, 1, Alloc_> base_type;
 26:     typedef typename base_type::point_type point_type;
 27:     typedef typename base_type::value_type value_type;
 28:     typedef typename base_type::atlas_type atlas_type;
 29:   protected:
 30:     int                       _localSize;
 31:     int                      *_offsets;
 32:     std::map<int, point_type> _invOrder;
 33:   public:
 34:     Numbering(MPI_Comm comm, const int debug = 0) : UniformSection<Point_, Value_, 1, Alloc_>(comm, debug), _localSize(0) {
 35:       this->_offsets    = new int[this->commSize()+1];
 36:       this->_offsets[0] = 0;
 37:     };
 38:     virtual ~Numbering() {
 39:       delete [] this->_offsets;
 40:     };
 41:   public: // Sizes
 42:     int        getLocalSize() const {return this->_localSize;};
 43:     void       setLocalSize(const int size) {this->_localSize = size;};
 44:     int        getGlobalSize() const {return this->_offsets[this->commSize()];};
 45:     int        getGlobalOffset(const int p) const {return this->_offsets[p];};
 46:     const int *getGlobalOffsets() const {return this->_offsets;};
 47:     void       setGlobalOffsets(const int offsets[]) {
 48:       for(int p = 0; p <= this->commSize(); ++p) {
 49:         this->_offsets[p] = offsets[p];
 50:       }
 51:     };
 52:   public: // Indices
 53:     virtual int getIndex(const point_type& point) {
 54:       const value_type& idx = this->restrictPoint(point)[0];
 55:       if (idx >= 0) {
 56:         return idx;
 57:       }
 58:       return -(idx+1);
 59:     };
 60:     virtual void setIndex(const point_type& point, const int index) {this->updatePoint(point, &index);};
 61:     virtual bool isLocal(const point_type& point) {return this->restrictPoint(point)[0] >= 0;};
 62:     virtual bool isRemote(const point_type& point) {return this->restrictPoint(point)[0] < 0;};
 63:     point_type getPoint(const int& index) {return this->_invOrder[index];};
 64:     void setPoint(const int& index, const point_type& point) {this->_invOrder[index] = point;};
 65:   };
 66:   template<typename Point_, typename Value_ = ALE::Point>
 67:   class GlobalOrder : public UniformSection<Point_, Value_> {
 68:   public:
 69:     typedef UniformSection<Point_, Value_> base_type;
 70:     typedef typename base_type::point_type point_type;
 71:     typedef typename base_type::value_type value_type;
 72:     typedef typename base_type::atlas_type atlas_type;
 73:   protected:
 74:     int  _localSize;
 75:     int *_offsets;
 76:   public:
 77:     GlobalOrder(MPI_Comm comm, const int debug = 0) : UniformSection<Point_, Value_>(comm, debug), _localSize(0) {
 78:       this->_offsets    = new int[this->commSize()+1];
 79:       this->_offsets[0] = 0;
 80:     };
 81:     ~GlobalOrder() {
 82:       delete [] this->_offsets;
 83:     };
 84:   public: // Sizes
 85:     int        getLocalSize() const {return this->_localSize;};
 86:     void       setLocalSize(const int size) {this->_localSize = size;};
 87:     int        getGlobalSize() const {return this->_offsets[this->commSize()];};
 88:     int        getGlobalOffset(const int p) const {return this->_offsets[p];};
 89:     const int *getGlobalOffsets() const {return this->_offsets;};
 90:     void       setGlobalOffsets(const int offsets[]) {
 91:       for(int p = 0; p <= this->commSize(); ++p) {
 92:         this->_offsets[p] = offsets[p];
 93:       }
 94:     };
 95:   public: // Indices
 96:     virtual int getIndex(const point_type& p) {
 97:       const int idx = this->restrictPoint(p)[0].prefix;
 98:       if (idx >= 0) {
 99:         return idx;
100:       }
101:       return -(idx+1);
102:     };
103:     virtual void setIndex(const point_type& p, const int index) {
104:       const value_type idx(index, this->restrictPoint(p)[0].index);
105:       this->updatePoint(p, &idx);
106:     };
107:     virtual bool isLocal(const point_type& p) {return this->restrictPoint(p)[0].prefix >= 0;};
108:     virtual bool isRemote(const point_type& p) {return this->restrictPoint(p)[0].prefix < 0;};
109:   };
110:   template<typename Bundle_, typename Value_ = int, typename Alloc_ = typename Bundle_::alloc_type>
111:   class NumberingFactory : ALE::ParallelObject {
112:   public:
113:     typedef Bundle_                                         bundle_type;
114:     typedef Alloc_                                          alloc_type;
115:     typedef Value_                                          value_type;
116:     typedef typename bundle_type::sieve_type                sieve_type;
117:     typedef typename bundle_type::point_type                point_type;
118:     typedef typename bundle_type::rank_type                 rank_type;
119:     typedef typename bundle_type::send_overlap_type         send_overlap_type;
120:     typedef typename bundle_type::recv_overlap_type         recv_overlap_type;
121:     typedef Numbering<point_type, value_type, alloc_type>   numbering_type;
122:     typedef typename alloc_type::template rebind<value_type>::other                              value_alloc_type;
123:     typedef std::map<bundle_type*, std::map<std::string, std::map<int, Obj<numbering_type> > > > numberings_type;
124:     typedef GlobalOrder<point_type>                         order_type;
125:     typedef typename order_type::value_type                 oValue_type;
126:     typedef typename alloc_type::template rebind<oValue_type>::other         oValue_alloc_type;
127:     typedef std::map<bundle_type*, std::map<std::string, Obj<order_type> > > orders_type;
128:   protected:
129:     numberings_type   _localNumberings;
130:     numberings_type   _numberings;
131:     orders_type       _orders;
132:     const value_type  _unknownNumber;
133:     const oValue_type _unknownOrder;
134:   protected:
135:     NumberingFactory(MPI_Comm comm, const int debug = 0) : ALE::ParallelObject(comm, debug), _unknownNumber(-1), _unknownOrder(-1, 0) {};
136:   public:
137:     ~NumberingFactory() {};
138:   public:
139:     static const Obj<NumberingFactory>& singleton(MPI_Comm comm, const int debug, bool cleanup = false) {
140:       static Obj<NumberingFactory> *_singleton = NULL;

142:       if (cleanup) {
143:         if (debug) {std::cout << "Destroying NumberingFactory" << std::endl;}
144:         if (_singleton) {delete _singleton;}
145:         _singleton = NULL;
146:       } else if (_singleton == NULL) {
147:         if (debug) {std::cout << "Creating new NumberingFactory" << std::endl;}
148:         _singleton  = new Obj<NumberingFactory>();
149:         *_singleton = new NumberingFactory(comm, debug);
150:       }
151:       return *_singleton;
152:     };
153:     void clear() {
154:       this->_localNumberings.clear();
155:       this->_numberings.clear();
156:       this->_orders.clear();
157:     };
158:   public: // Dof ordering
159:     template<typename Section_>
160:     void orderPointNew(const Obj<Section_>& section, const Obj<sieve_type>& sieve, const typename Section_::point_type& point, value_type& offset, value_type& bcOffset, const Obj<send_overlap_type>& sendOverlap = NULL) {
161:       const typename Section_::chart_type& chart = section->getChart();
162:       int&                                 idx   = section->getIndex(point);

164:       // If the point does not exist in the chart, throw an error
165:       if (chart.count(point) == 0) {
166:         throw ALE::Exception("Unknown point in ordering");
167:       }
168:       // If the point has not been ordered
169:       if (idx == -1) {
170:         // Recurse to its cover
171:         const Obj<typename sieve_type::coneSequence>& cone = sieve->cone(point);
172:         typename sieve_type::coneSequence::iterator   end  = cone->end();

174:         for(typename sieve_type::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
175:           if (this->_debug > 1) {std::cout << "    Recursing to " << *c_iter << std::endl;}
176:           this->orderPoint(section, sieve, *c_iter, offset, bcOffset, sendOverlap);
177:         }
178:         const int dim  = section->getFiberDimension(point);
179:         const int cDim = section->getConstraintDimension(point);
180:         const int fDim = dim - cDim;

182:         // If the point has constrained variables
183:         if (cDim) {
184:           if (this->_debug > 1) {std::cout << "  Ordering boundary point " << point << " at " << bcOffset << std::endl;}
185:           section->setIndexBC(point, bcOffset);
186:           bcOffset += cDim;
187:         }
188:         // If the point has free variables
189:         if (fDim) {
190:           bool number = true;

192:           // Maybe use template specialization here
193:           if (!sendOverlap.isNull() && sendOverlap->capContains(point)) {
194:             const Obj<typename send_overlap_type::supportSequence>& ranks = sendOverlap->support(point);

196:             for(typename send_overlap_type::supportSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
197:               if (this->commRank() > *r_iter) {
198:                 number = false;
199:                 break;
200:               }
201:             }
202:           }
203:           if (number) {
204:             if (this->_debug > 1) {std::cout << "  Ordering point " << point << " at " << offset << std::endl;}
205:             section->setIndex(point, offset);
206:             offset += dim;
207:           } else {
208:             if (this->_debug > 1) {std::cout << "  Ignoring ghost point " << point << std::endl;}
209:           }
210:         }
211:       }
212:     };
213:     template<typename Section_>
214:     void orderPoint(const Obj<Section_>& section, const Obj<sieve_type>& sieve, const typename Section_::point_type& point, value_type& offset, value_type& bcOffset, const Obj<send_overlap_type>& sendOverlap = NULL) {
215:       const Obj<typename Section_::atlas_type>&     atlas = section->getAtlas();
216:       const Obj<typename sieve_type::coneSequence>& cone = sieve->cone(point);
217:       typename sieve_type::coneSequence::iterator   end  = cone->end();
218:       typename Section_::index_type                 idx  = section->getAtlas()->restrictPoint(point)[0];
219:       const value_type&                             dim  = idx.prefix;
220:       const typename Section_::index_type           defaultIdx(0, -1);

222:       if (atlas->getChart().count(point) == 0) {
223:         idx = defaultIdx;
224:       }
225:       if (idx.index == -1) {
226:         for(typename sieve_type::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
227:           if (this->_debug > 1) {std::cout << "    Recursing to " << *c_iter << std::endl;}
228:           this->orderPoint(section, sieve, *c_iter, offset, bcOffset, sendOverlap);
229:         }
230:         if (dim > 0) {
231:           bool number = true;

233:           // Maybe use template specialization here
234:           if (!sendOverlap.isNull() && sendOverlap->capContains(point)) {
235:             const Obj<typename send_overlap_type::supportSequence>& ranks = sendOverlap->support(point);

237:             for(typename send_overlap_type::supportSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
238:               if (this->commRank() > *r_iter) {
239:                 number = false;
240:                 break;
241:               }
242:             }
243:           }
244:           if (number) {
245:             if (this->_debug > 1) {std::cout << "  Ordering point " << point << " at " << offset << std::endl;}
246:             idx.index = offset;
247:             atlas->updatePoint(point, &idx);
248:             offset += dim;
249:           } else {
250:             if (this->_debug > 1) {std::cout << "  Ignoring ghost point " << point << std::endl;}
251:           }
252:         } else if (dim < 0) {
253:           if (this->_debug > 1) {std::cout << "  Ordering boundary point " << point << " at " << bcOffset << std::endl;}
254:           idx.index = bcOffset;
255:           atlas->updatePoint(point, &idx);
256:           bcOffset += dim;
257:         }
258:       }
259:     };
260:     template<typename Section_>
261:     void orderPatch(const Obj<Section_>& section, const Obj<sieve_type>& sieve, const Obj<send_overlap_type>& sendOverlap = NULL, const value_type offset = 0, const value_type bcOffset = -2) {
262:       const typename Section_::chart_type& chart = section->getChart();
263:       int off   = offset;
264:       int bcOff = bcOffset;

266:       if (this->_debug > 1) {std::cout << "Ordering patch" << std::endl;}
267:       for(typename Section_::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
268:         if (this->_debug > 1) {std::cout << "Ordering closure of point " << *p_iter << std::endl;}
269:         this->orderPoint(section, sieve, *p_iter, off, bcOff, sendOverlap);
270:       }
271:       for(typename Section_::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
272:         const int& idx  = section->getIndex(*p_iter);

274:         if (idx < 0) {
275:           if (this->_debug > 1) {std::cout << "Correcting boundary offset of point " << *p_iter << std::endl;}
276:           section->setIndex(*p_iter, off - (idx + 2));
277:         }
278:       }
279:     };
280:   public: // Numbering
281:     // Number all local points
282:     //   points in the overlap are only numbered by the owner with the lowest rank
283:     template<typename Sequence_>
284:     void constructLocalNumbering(const Obj<numbering_type>& numbering, const Obj<send_overlap_type>& sendOverlap, const Obj<Sequence_>& points) {
285:       const int debug = sendOverlap->debug();
286:       int localSize = 0;

288:       if (debug) {std::cout << "["<<numbering->commRank()<<"] Constructing local numbering" << std::endl;}
289:       numbering->setFiberDimension(points, 1);
290:       for(typename Sequence_::iterator l_iter = points->begin(); l_iter != points->end(); ++l_iter) {
291:         value_type val;

293:         if (debug) {std::cout << "["<<numbering->commRank()<<"]   Checking point " << *l_iter << std::endl;}
294:         if (sendOverlap->capContains(*l_iter)) {
295:           const Obj<typename send_overlap_type::traits::supportSequence>& sendPatches = sendOverlap->support(*l_iter);
296:           int minRank = sendOverlap->commSize();

298:           for(typename send_overlap_type::traits::supportSequence::iterator p_iter = sendPatches->begin(); p_iter != sendPatches->end(); ++p_iter) {
299:             if (*p_iter < minRank) minRank = *p_iter;
300:           }
301:           if (minRank < sendOverlap->commRank()) {
302:             if (debug) {std::cout << "["<<numbering->commRank()<<"]     remote point, on proc " << minRank << std::endl;}
303:             val = this->_unknownNumber;
304:           } else {
305:             if (debug) {std::cout << "["<<numbering->commRank()<<"]     local point" << std::endl;}
306:             val = localSize++;
307:           }
308:         } else {
309:           if (debug) {std::cout << "["<<numbering->commRank()<<"]     local point" << std::endl;}
310:           val = localSize++;
311:         }
312:         if (debug) {std::cout << "["<<numbering->commRank()<<"]     has number " << val << std::endl;}
313:         numbering->updatePoint(*l_iter, &val);
314:       }
315:       if (debug) {std::cout << "["<<numbering->commRank()<<"]   local points" << std::endl;}
316:       numbering->setLocalSize(localSize);
317:     };
318:     // Order all local points
319:     //   points in the overlap are only ordered by the owner with the lowest rank
320:     template<typename Sequence_, typename Section_>
321:     void constructLocalOrder(const Obj<order_type>& order, const Obj<send_overlap_type>& sendOverlap, const Sequence_& points, const Obj<Section_>& section) {
322:       const int debug = sendOverlap->debug();
323:       int localSize = 0;

325:       if (debug) {std::cout << "["<<order->commRank()<<"] Constructing local ordering" << std::endl;}
326:       for(typename Sequence_::const_iterator l_iter = points.begin(); l_iter != points.end(); ++l_iter) {
327:         order->setFiberDimension(*l_iter, 1);
328:       }
329:       for(typename Sequence_::const_iterator l_iter = points.begin(); l_iter != points.end(); ++l_iter) {
330:         oValue_type val;

332:         if (debug) {std::cout << "["<<order->commRank()<<"]   Checking point " << *l_iter << std::endl;}
333:         if (sendOverlap->capContains(*l_iter)) {
334:           const Obj<typename send_overlap_type::traits::supportSequence>& sendPatches = sendOverlap->support(*l_iter);
335:           int minRank = sendOverlap->commSize();

337:           for(typename send_overlap_type::traits::supportSequence::iterator p_iter = sendPatches->begin(); p_iter != sendPatches->end(); ++p_iter) {
338:             if (*p_iter < minRank) minRank = *p_iter;
339:           }
340:           if (minRank < sendOverlap->commRank()) {
341:             if (debug) {std::cout << "["<<order->commRank()<<"]     remote point, on proc " << minRank << std::endl;}
342:             val = this->_unknownOrder;
343:           } else {
344:             if (debug) {std::cout << "["<<order->commRank()<<"]     local point" << std::endl;}
345:             val.prefix = localSize;
346:             val.index  = section->getConstrainedFiberDimension(*l_iter);
347:           }
348:         } else {
349:           if (debug) {std::cout << "["<<order->commRank()<<"]     local point" << std::endl;}
350:           val.prefix = localSize;
351:           val.index  = section->getConstrainedFiberDimension(*l_iter);
352:         }
353:         if (debug) {std::cout << "["<<order->commRank()<<"]     has offset " << val.prefix << " and size " << val.index << std::endl;}
354:         localSize += val.index;
355:         order->updatePoint(*l_iter, &val);
356:       }
357:       if (debug) {std::cout << "["<<order->commRank()<<"]   local size" << localSize << std::endl;}
358:       order->setLocalSize(localSize);
359:     };
360:     // Calculate process offsets
361:     template<typename Numbering>
362:     void calculateOffsets(const Obj<Numbering>& numbering) {
363:       int  localSize = numbering->getLocalSize();
364:       int *offsets   = new int[numbering->commSize()+1];

366:       offsets[0] = 0;
367:       MPI_Allgather(&localSize, 1, MPI_INT, &(offsets[1]), 1, MPI_INT, numbering->comm());
368:       for(int p = 2; p <= numbering->commSize(); p++) {
369:         offsets[p] += offsets[p-1];
370:       }
371:       numbering->setGlobalOffsets(offsets);
372:       delete [] offsets;
373:     };
374:     // Update local offsets based upon process offsets
375:     //   TODO: The sequence should be const, but LabelSifter has no proper const_iterator
376:     template<typename Numbering, typename Sequence>
377:     void updateOrder(const Obj<Numbering>& numbering, Sequence& points) {
378:       const typename Numbering::value_type val = numbering->getGlobalOffset(numbering->commRank());

380:       for(typename Sequence::const_iterator l_iter = points.begin(); l_iter != points.end(); ++l_iter) {
381:         if (numbering->isLocal(*l_iter)) {
382:           numbering->updateAddPoint(*l_iter, &val);
383:         }
384:       }
385:     };
386:     // Communicate numbers in the overlap
387:     void completeNumbering(const Obj<numbering_type>& numbering, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, bool allowDuplicates = false) {
388:       typedef Field<send_overlap_type, int, Section<point_type, value_type, value_alloc_type> > send_section_type;
389:       typedef Field<recv_overlap_type, int, Section<point_type, value_type, value_alloc_type> > recv_section_type;
390:       typedef typename ALE::DiscreteSieve<point_type, alloc_type>                   dsieve_type;
391:       typedef typename ALE::Topology<int, dsieve_type, alloc_type>                  dtopology_type;
392:       typedef typename ALE::New::SectionCompletion<dtopology_type, int, alloc_type> completion;
393:       const Obj<send_section_type> sendSection = new send_section_type(numbering->comm(), this->debug());
394:       const Obj<recv_section_type> recvSection = new recv_section_type(numbering->comm(), sendSection->getTag(), this->debug());
395:       const int debug = sendOverlap->debug();

397:       if (debug) {std::cout << "["<<numbering->commRank()<<"] Completing numbering" << std::endl;}
398:       completion::completeSection(sendOverlap, recvOverlap, numbering->getAtlas(), numbering, sendSection, recvSection);
399: #if 1
400:       const Obj<typename recv_overlap_type::traits::baseSequence> rPoints = recvOverlap->base();

402:       for(typename recv_overlap_type::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rPoints->end(); ++p_iter) {
403:         const Obj<typename recv_overlap_type::coneSequence>& ranks      = recvOverlap->cone(*p_iter);
404:         const typename recv_overlap_type::target_type&       localPoint = *p_iter;

406:         for(typename recv_overlap_type::coneSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
407:           const typename recv_overlap_type::target_type&       remotePoint = r_iter.color();
408:           const int                                            rank        = *r_iter;
409:           const Obj<typename recv_section_type::section_type>& section     = recvSection->getSection(rank);
410:           const typename recv_section_type::value_type        *values      = section->restrictPoint(remotePoint);

412:           if (section->getFiberDimension(remotePoint) == 0) continue;
413:           if (debug) {std::cout << "["<<numbering->commRank()<<"]     local point " << localPoint << " remote point " << remotePoint << " number " << values[0] << std::endl;}
414:           if (values[0] >= 0) {
415:             if (debug) {std::cout << "["<<numbering->commRank()<<"] local point " << localPoint << " dim " << numbering->getAtlas()->getFiberDimension(localPoint) << std::endl;}
416:             if (numbering->isLocal(localPoint) && !allowDuplicates) {
417:               ostringstream msg;
418:               msg << "["<<numbering->commRank()<<"]Multiple indices for local point " << localPoint << " remote point " << remotePoint << " from " << rank << " with index " << values[0];
419:               throw ALE::Exception(msg.str().c_str());
420:             }
421:             if (numbering->getAtlas()->getFiberDimension(localPoint) == 0) {
422:               ostringstream msg;
423:               msg << "["<<numbering->commRank()<<"]Unexpected local point " << localPoint << " remote point " << remotePoint << " from " << rank << " with index " << values[0];
424:               throw ALE::Exception(msg.str().c_str());
425:             }
426:             int val = -(values[0]+1);
427:             numbering->updatePoint(localPoint, &val);
428:           }
429:         }
430:       }
431: #else
432:       const typename recv_section_type::sheaf_type& patches = recvSection->getPatches();

434:       for(typename recv_section_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
435:         const typename recv_section_type::patch_type&        rPatch  = p_iter->first;
436:         const Obj<typename recv_section_type::section_type>& section = recvSection->getSection(rPatch);
437:         const typename recv_section_type::chart_type&        points  = section->getChart();

439:         for(typename recv_section_type::chart_type::const_iterator r_iter = points.begin(); r_iter != points.end(); ++r_iter) {
440:           const typename recv_section_type::point_type& point  = *r_iter;
441:           const typename recv_section_type::value_type *values = section->restrictPoint(point);

443:           if (section->getFiberDimension(point) == 0) continue;
444:           if (values[0] >= 0) {
445:             if (numbering->isLocal(point) && !allowDuplicates) {
446:               ostringstream msg;
447:               msg << "["<<numbering->commRank()<<"]Multiple indices for point " << point << " from " << rPatch << " with index " << values[0];
448:               throw ALE::Exception(msg.str().c_str());
449:             }
450:             if (numbering->getAtlas()->getFiberDimension(point) == 0) {
451:               ostringstream msg;
452:               msg << "["<<numbering->commRank()<<"]Unexpected point " << point << " from " << rPatch << " with index " << values[0];
453:               throw ALE::Exception(msg.str().c_str());
454:             }
455:             int val = -(values[0]+1);
456:             numbering->updatePoint(point, &val);
457:           }
458:         }
459:       }
460: #endif
461:     };
462:     // Communicate (size,offset)s in the overlap
463:     void completeOrder(const Obj<order_type>& order, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, bool allowDuplicates = false) {
464:       typedef Field<send_overlap_type, int, Section<point_type, oValue_type, oValue_alloc_type> > send_section_type;
465:       typedef Field<recv_overlap_type, int, Section<point_type, oValue_type, oValue_alloc_type> > recv_section_type;
466:       typedef ConstantSection<point_type, int, alloc_type>                          constant_sizer;
467:       typedef typename ALE::DiscreteSieve<point_type, alloc_type>                   dsieve_type;
468:       typedef typename ALE::Topology<int, dsieve_type, alloc_type>                  dtopology_type;
469:       typedef typename ALE::New::SectionCompletion<dtopology_type, int, alloc_type> completion;
470:       const Obj<send_section_type> sendSection = new send_section_type(order->comm(), this->debug());
471:       const Obj<recv_section_type> recvSection = new recv_section_type(order->comm(), sendSection->getTag(), this->debug());
472:       const int debug = sendOverlap->debug();

474:       if (debug) {std::cout << "["<<order->commRank()<<"] Completing ordering" << std::endl;}
475:       completion::completeSection(sendOverlap, recvOverlap, order->getAtlas(), order, sendSection, recvSection);
476:       Obj<typename recv_overlap_type::traits::baseSequence> recvPoints = recvOverlap->base();

478:       for(typename recv_overlap_type::traits::baseSequence::iterator p_iter = recvPoints->begin(); p_iter != recvPoints->end(); ++p_iter) {
479:         if (!order->hasPoint(*p_iter)) {
480:           order->setFiberDimension(*p_iter, 1);
481:           order->updatePoint(*p_iter, &this->_unknownOrder);
482:         }
483:       }
484:       for(typename recv_overlap_type::traits::baseSequence::iterator p_iter = recvPoints->begin(); p_iter != recvPoints->end(); ++p_iter) {
485:         const Obj<typename recv_overlap_type::traits::coneSequence>& ranks      = recvOverlap->cone(*p_iter);
486:         const typename recv_overlap_type::target_type&               localPoint = *p_iter;

488:         for(typename recv_overlap_type::traits::coneSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
489:           const typename recv_overlap_type::target_type&       remotePoint = r_iter.color();
490:           const int                                            rank        = *r_iter;
491:           const Obj<typename recv_section_type::section_type>& section     = recvSection->getSection(rank);
492:           const typename recv_section_type::value_type        *values      = section->restrictPoint(remotePoint);

494:           if (section->getFiberDimension(remotePoint) == 0) continue;
495:           if (debug) {std::cout << "["<<order->commRank()<<"]     local point " << localPoint << " remote point " << remotePoint<<"("<<rank<<")" << " offset " << values[0].prefix << " and size " << values[0].index << std::endl;}
496:           if (values[0].index == 0) continue;
497:           if (values[0].prefix >= 0) {
498:             if (order->isLocal(localPoint)) {
499:               if (!allowDuplicates) {
500:                 ostringstream msg;
501:                 msg << "["<<order->commRank()<<"]Multiple indices for local point " << localPoint << " remote point " << remotePoint << " from " << rank << " with index " << values[0];
502:                 throw ALE::Exception(msg.str().c_str());
503:               }
504:               continue;
505:             }
506:             const oValue_type val(-(values[0].prefix+1), values[0].index);
507:             order->updatePoint(localPoint, &val);
508:           } else {
509:             if (order->isLocal(localPoint)) continue;
510:             order->updatePoint(localPoint, values);
511:           }
512:         }
513:       }
514:     };
515:     // Construct a full global numbering
516:     template<typename Sequence>
517:     void constructNumbering(const Obj<numbering_type>& numbering, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Obj<Sequence>& points) {
518:       this->constructLocalNumbering(numbering, sendOverlap, points);
519:       this->calculateOffsets(numbering);
520:       this->updateOrder(numbering, *points.ptr());
521:       this->completeNumbering(numbering, sendOverlap, recvOverlap);
522:     };
523:     // Construct a full global order
524:     template<typename Sequence, typename Section>
525:     void constructOrder(const Obj<order_type>& order, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Sequence& points, const Obj<Section>& section) {
526:       this->constructLocalOrder(order, sendOverlap, points, section);
527:       this->calculateOffsets(order);
528:       this->updateOrder(order, points);
529:       this->completeOrder(order, sendOverlap, recvOverlap);
530:     };
531:     template<typename Sequence, typename Section>
532:     void constructOrder(const Obj<order_type>& order, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Obj<Sequence>& points, const Obj<Section>& section) {
533:       this->constructLocalOrder(order, sendOverlap, *points.ptr(), section);
534:       this->calculateOffsets(order);
535:       this->updateOrder(order, *points.ptr());
536:       this->completeOrder(order, sendOverlap, recvOverlap);
537:     };
538:   public:
539:     // Construct the inverse map from numbers to points
540:     //   If we really need this, then we should consider using a label
541:     void constructInverseOrder(const Obj<numbering_type>& numbering) {
542:       const typename numbering_type::chart_type& chart = numbering->getChart();

544:       for(typename numbering_type::chart_type::iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
545:         numbering->setPoint(numbering->getIndex(*p_iter), *p_iter);
546:       }
547:     };
548:   public: // Real interface
549:     template<typename ABundle_>
550:     const Obj<numbering_type>& getLocalNumbering(const Obj<ABundle_>& bundle, const int depth) {
551:       if ((this->_localNumberings.find(bundle.ptr()) == this->_localNumberings.end()) ||
552:           (this->_localNumberings[bundle.ptr()].find("depth") == this->_localNumberings[bundle.ptr()].end()) ||
553:           (this->_localNumberings[bundle.ptr()]["depth"].find(depth) == this->_localNumberings[bundle.ptr()]["depth"].end())) {
554:         Obj<numbering_type>    numbering   = new numbering_type(bundle->comm(), bundle->debug());
555:         Obj<send_overlap_type> sendOverlap = new send_overlap_type(bundle->comm(), bundle->debug());

557:         this->constructLocalNumbering(numbering, sendOverlap, bundle->depthStratum(depth));
558:         if (this->_debug) {std::cout << "Creating new local numbering: ptr " << bundle.ptr() << " depth " << depth << std::endl;}
559:         this->_localNumberings[bundle.ptr()]["depth"][depth] = numbering;
560:       } else {
561:         if (this->_debug) {std::cout << "Using old local numbering: ptr " << bundle.ptr() << " depth " << depth << std::endl;}
562:       }
563:       return this->_localNumberings[bundle.ptr()]["depth"][depth];
564:     };
565:     template<typename ABundle_>
566:     const Obj<numbering_type>& getNumbering(const Obj<ABundle_>& bundle, const int depth) {
567:       if ((this->_numberings.find(bundle.ptr()) == this->_numberings.end()) ||
568:           (this->_numberings[bundle.ptr()].find("depth") == this->_numberings[bundle.ptr()].end()) ||
569:           (this->_numberings[bundle.ptr()]["depth"].find(depth) == this->_numberings[bundle.ptr()]["depth"].end())) {
570:         bundle->constructOverlap();
571:         Obj<numbering_type>    numbering   = new numbering_type(bundle->comm(), bundle->debug());
572:         Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
573:         Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();

575: //         std::cout << "["<<bundle->commRank()<<"]Creating new numbering: fixed depth value " << depth << std::endl;
576:         this->constructNumbering(numbering, sendOverlap, recvOverlap, bundle->depthStratum(depth));
577:         if (this->_debug) {std::cout << "Creating new numbering: depth " << depth << std::endl;}
578:         this->_numberings[bundle.ptr()]["depth"][depth] = numbering;
579: //       } else {
580: //         std::cout << "["<<bundle->commRank()<<"]Using old numbering: fixed depth value " << depth << std::endl;
581:       }
582:       return this->_numberings[bundle.ptr()]["depth"][depth];
583:     };
584:     template<typename ABundle_>
585:     const Obj<numbering_type>& getNumbering(const Obj<ABundle_>& bundle, const std::string& labelname, const int value) {
586:       if ((this->_numberings.find(bundle.ptr()) == this->_numberings.end()) ||
587:           (this->_numberings[bundle.ptr()].find(labelname) == this->_numberings[bundle.ptr()].end()) ||
588:           (this->_numberings[bundle.ptr()][labelname].find(value) == this->_numberings[bundle.ptr()][labelname].end())) {
589:         bundle->constructOverlap();
590:         Obj<numbering_type>    numbering   = new numbering_type(bundle->comm(), bundle->debug());
591:         Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
592:         Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();

594:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Creating new numbering: " << labelname << " value " << value << std::endl;}
595:         this->constructNumbering(numbering, sendOverlap, recvOverlap, bundle->getLabelStratum(labelname, value));
596:         this->_numberings[bundle.ptr()][labelname][value] = numbering;
597:       } else {
598:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Using old numbering: " << labelname << " value " << value << std::endl;}
599:       }
600:       return this->_numberings[bundle.ptr()][labelname][value];
601:     };
602:     template<typename ABundle_, typename Section_>
603:     const Obj<order_type>& getGlobalOrder(const Obj<ABundle_>& bundle, const std::string& name, const Obj<Section_>& section) {
604:       if ((this->_orders.find(bundle.ptr()) == this->_orders.end()) ||
605:           (this->_orders[bundle.ptr()].find(name) == this->_orders[bundle.ptr()].end())) {
606:         bundle->constructOverlap();
607:         Obj<order_type>        order       = new order_type(bundle->comm(), bundle->debug());
608:         Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
609:         Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();

611:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Creating new global order: " << name << std::endl;}
612:         this->constructOrder(order, sendOverlap, recvOverlap, section->getChart(), section);
613:         this->_orders[bundle.ptr()][name] = order;
614:       } else {
615:         if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Using old global order: " << name << std::endl;}
616:       }
617:       return this->_orders[bundle.ptr()][name];
618:     };
619:     template<typename ABundle_>
620:     void setGlobalOrder(const Obj<ABundle_>& bundle, const std::string& name, const Obj<order_type>& order) {
621:       this->_orders[bundle.ptr()][name] = order;
622:     };
623:   };
624: }
625: #endif