Actual source code: Numbering.hh
petsc-3.3-p7 2013-05-11
1: #ifndef included_ALE_Numbering_hh
2: #define included_ALE_Numbering_hh
4: #ifndef included_ALE_ParallelMapping_hh
5: #include <sieve/ParallelMapping.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 : public 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::label_type label_type;
118: typedef typename bundle_type::point_type point_type;
119: typedef typename bundle_type::rank_type rank_type;
120: typedef typename bundle_type::send_overlap_type send_overlap_type;
121: typedef typename bundle_type::recv_overlap_type recv_overlap_type;
122: typedef Numbering<point_type, value_type, alloc_type> numbering_type;
123: typedef typename alloc_type::template rebind<value_type>::other value_alloc_type;
124: typedef std::map<bundle_type*, std::map<std::string, std::map<int, Obj<numbering_type> > > > numberings_type;
125: typedef GlobalOrder<point_type> order_type;
126: typedef typename order_type::value_type oValue_type;
127: typedef typename alloc_type::template rebind<oValue_type>::other oValue_alloc_type;
128: typedef std::map<bundle_type*, std::map<std::string, Obj<order_type> > > orders_type;
129: protected:
130: numberings_type _localNumberings;
131: numberings_type _numberings;
132: orders_type _localOrders;
133: orders_type _orders;
134: orders_type _ordersBC;
135: const value_type _unknownNumber;
136: const oValue_type _unknownOrder;
137: protected:
138: NumberingFactory(MPI_Comm comm, const int debug = 0) : ALE::ParallelObject(comm, debug), _unknownNumber(-1), _unknownOrder(-1, 0) {};
139: public:
140: ~NumberingFactory() {};
141: public:
142: static const Obj<NumberingFactory>& singleton(MPI_Comm comm, const int debug, bool cleanup = false) {
143: static Obj<NumberingFactory> *_singleton = NULL;
145: if (cleanup) {
146: if (debug) {std::cout << "Destroying NumberingFactory" << std::endl;}
147: if (_singleton) {delete _singleton;}
148: _singleton = NULL;
149: } else if (_singleton == NULL) {
150: if (debug) {std::cout << "Creating new NumberingFactory" << std::endl;}
151: _singleton = new Obj<NumberingFactory>();
152: *_singleton = new NumberingFactory(comm, debug);
153: }
154: return *_singleton;
155: };
156: void clear() {
157: this->_localNumberings.clear();
158: this->_numberings.clear();
159: this->_localOrders.clear();
160: this->_orders.clear();
161: this->_ordersBC.clear();
162: };
163: public: // Dof ordering
164: template<typename Section_>
165: 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) {
166: const typename Section_::chart_type& chart = section->getChart();
167: int& idx = section->getIndex(point);
169: // If the point does not exist in the chart, throw an error
170: if (chart.count(point) == 0) {
171: throw ALE::Exception("Unknown point in ordering");
172: }
173: // If the point has not been ordered
174: if (idx == -1) {
175: // Recurse to its cover
176: const Obj<typename sieve_type::coneSequence>& cone = sieve->cone(point);
177: typename sieve_type::coneSequence::iterator end = cone->end();
179: for(typename sieve_type::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
180: if (this->_debug > 1) {std::cout << " Recursing to " << *c_iter << std::endl;}
181: this->orderPoint(section, sieve, *c_iter, offset, bcOffset, sendOverlap);
182: }
183: const int dim = section->getFiberDimension(point);
184: const int cDim = section->getConstraintDimension(point);
185: const int fDim = dim - cDim;
187: // If the point has constrained variables
188: if (cDim) {
189: if (this->_debug > 1) {std::cout << " Ordering boundary point " << point << " at " << bcOffset << std::endl;}
190: section->setIndexBC(point, bcOffset);
191: bcOffset += cDim;
192: }
193: // If the point has free variables
194: if (fDim) {
195: bool number = true;
197: // Maybe use template specialization here
198: if (!sendOverlap.isNull() && sendOverlap->capContains(point)) {
199: const Obj<typename send_overlap_type::supportSequence>& ranks = sendOverlap->support(point);
201: for(typename send_overlap_type::supportSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
202: if (this->commRank() > *r_iter) {
203: number = false;
204: break;
205: }
206: }
207: }
208: if (number) {
209: if (this->_debug > 1) {std::cout << " Ordering point " << point << " at " << offset << std::endl;}
210: section->setIndex(point, offset);
211: offset += dim;
212: } else {
213: if (this->_debug > 1) {std::cout << " Ignoring ghost point " << point << std::endl;}
214: }
215: }
216: }
217: }
218: template<typename Section_>
219: 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) {
220: const Obj<typename Section_::atlas_type>& atlas = section->getAtlas();
221: const Obj<typename sieve_type::coneSequence>& cone = sieve->cone(point);
222: typename sieve_type::coneSequence::iterator end = cone->end();
223: typename Section_::index_type idx = section->getAtlas()->restrictPoint(point)[0];
224: const value_type& dim = idx.prefix;
225: const typename Section_::index_type defaultIdx(0, -1);
227: if (atlas->getChart().count(point) == 0) {
228: idx = defaultIdx;
229: }
230: if (idx.index == -1) {
231: for(typename sieve_type::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
232: if (this->_debug > 1) {std::cout << " Recursing to " << *c_iter << std::endl;}
233: this->orderPoint(section, sieve, *c_iter, offset, bcOffset, sendOverlap);
234: }
235: if (dim > 0) {
236: bool number = true;
238: // Maybe use template specialization here
239: if (!sendOverlap.isNull() && sendOverlap->capContains(point)) {
240: const Obj<typename send_overlap_type::supportSequence>& ranks = sendOverlap->support(point);
242: for(typename send_overlap_type::supportSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
243: if (this->commRank() > *r_iter) {
244: number = false;
245: break;
246: }
247: }
248: }
249: if (number) {
250: if (this->_debug > 1) {std::cout << " Ordering point " << point << " at " << offset << std::endl;}
251: idx.index = offset;
252: atlas->updatePoint(point, &idx);
253: offset += dim;
254: } else {
255: if (this->_debug > 1) {std::cout << " Ignoring ghost point " << point << std::endl;}
256: }
257: } else if (dim < 0) {
258: if (this->_debug > 1) {std::cout << " Ordering boundary point " << point << " at " << bcOffset << std::endl;}
259: idx.index = bcOffset;
260: atlas->updatePoint(point, &idx);
261: bcOffset += dim;
262: }
263: }
264: }
265: template<typename Section_>
266: 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) {
267: const typename Section_::chart_type& chart = section->getChart();
268: int off = offset;
269: int bcOff = bcOffset;
271: if (this->_debug > 1) {std::cout << "Ordering patch" << std::endl;}
272: for(typename Section_::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
273: if (this->_debug > 1) {std::cout << "Ordering closure of point " << *p_iter << std::endl;}
274: this->orderPoint(section, sieve, *p_iter, off, bcOff, sendOverlap);
275: }
276: for(typename Section_::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
277: const int& idx = section->getIndex(*p_iter);
279: if (idx < 0) {
280: if (this->_debug > 1) {std::cout << "Correcting boundary offset of point " << *p_iter << std::endl;}
281: section->setIndex(*p_iter, off - (idx + 2));
282: }
283: }
284: }
285: public: // Numbering
286: // Number all local points
287: // points in the overlap are only numbered by the owner with the lowest rank
288: template<typename Iterator_>
289: void constructLocalNumbering(const Obj<numbering_type>& numbering, const Obj<send_overlap_type>& sendOverlap, const Iterator_& pointsBegin, const Iterator_& pointsEnd) {
290: const int debug = sendOverlap->debug();
291: int localSize = 0;
293: if (debug) {std::cout << "["<<numbering->commRank()<<"] Constructing local numbering" << std::endl;}
294: for(Iterator_ l_iter = pointsBegin; l_iter != pointsEnd; ++l_iter) {
295: numbering->setFiberDimension(*l_iter, 1);
296: }
297: for(Iterator_ l_iter = pointsBegin; l_iter != pointsEnd; ++l_iter) {
298: value_type val;
300: if (debug) {std::cout << "["<<numbering->commRank()<<"] Checking point " << *l_iter << std::endl;}
301: if (sendOverlap->capContains(*l_iter)) {
302: const Obj<typename send_overlap_type::supportSequence>& sendRanks = sendOverlap->support(*l_iter);
303: int minRank = sendOverlap->commSize();
305: for(typename send_overlap_type::supportSequence::iterator p_iter = sendRanks->begin(); p_iter != sendRanks->end(); ++p_iter) {
306: if (*p_iter < minRank) minRank = *p_iter;
307: }
308: if (minRank < sendOverlap->commRank()) {
309: if (debug) {std::cout << "["<<numbering->commRank()<<"] remote point, on proc " << minRank << std::endl;}
310: val = this->_unknownNumber;
311: } else {
312: if (debug) {std::cout << "["<<numbering->commRank()<<"] local point" << std::endl;}
313: val = localSize++;
314: }
315: } else {
316: if (debug) {std::cout << "["<<numbering->commRank()<<"] local point" << std::endl;}
317: val = localSize++;
318: }
319: if (debug) {std::cout << "["<<numbering->commRank()<<"] has number " << val << std::endl;}
320: numbering->updatePoint(*l_iter, &val);
321: }
322: if (debug) {std::cout << "["<<numbering->commRank()<<"] local points" << std::endl;}
323: numbering->setLocalSize(localSize);
324: }
325: // Order all local points
326: // points in the overlap are only ordered by the owner with the lowest rank
327: template<typename Iterator_, typename Section_>
328: void constructLocalOrder(const Obj<order_type>& order, const Obj<send_overlap_type>& sendOverlap, const Iterator_& pointsBegin, const Iterator_& pointsEnd, const Obj<Section_>& section, const int space = -1, const bool withBC = false, const Obj<label_type>& label = PETSC_NULL) {
329: const int debug = sendOverlap->debug();
330: int localSize = 0;
332: if (debug) {std::cout << "["<<order->commRank()<<"] Constructing local ordering" << std::endl;}
333: for(Iterator_ l_iter = pointsBegin; l_iter != pointsEnd; ++l_iter) {
334: order->setFiberDimension(*l_iter, 1);
335: }
336: for(Iterator_ l_iter = pointsBegin; l_iter != pointsEnd; ++l_iter) {
337: oValue_type val;
339: if (debug) {std::cout << "["<<order->commRank()<<"] Checking point " << *l_iter << std::endl;}
340: if (sendOverlap->capContains(*l_iter)) {
341: const Obj<typename send_overlap_type::supportSequence>& sendPatches = sendOverlap->support(*l_iter);
342: int minRank = sendOverlap->commSize();
344: for(typename send_overlap_type::supportSequence::iterator p_iter = sendPatches->begin(); p_iter != sendPatches->end(); ++p_iter) {
345: if (*p_iter < minRank) minRank = *p_iter;
346: }
347: bool remotePoint = (minRank < sendOverlap->commRank()) || (!label.isNull() && (label->cone(*l_iter)->size() > 0));
349: if (remotePoint) {
350: if (debug) {std::cout << "["<<order->commRank()<<"] remote point, on proc " << minRank << std::endl;}
351: val = this->_unknownOrder;
352: } else {
353: if (debug) {std::cout << "["<<order->commRank()<<"] local point" << std::endl;}
354: val.prefix = localSize;
355: if (withBC) {
356: val.index = space < 0 ? section->getFiberDimension(*l_iter) : section->getFiberDimension(*l_iter, space);
357: } else {
358: val.index = space < 0 ? section->getConstrainedFiberDimension(*l_iter) : section->getConstrainedFiberDimension(*l_iter, space);
359: }
360: }
361: } else {
362: if (debug) {std::cout << "["<<order->commRank()<<"] local point" << std::endl;}
363: val.prefix = localSize;
364: if (withBC) {
365: val.index = space < 0 ? section->getFiberDimension(*l_iter) : section->getFiberDimension(*l_iter, space);
366: } else {
367: val.index = space < 0 ? section->getConstrainedFiberDimension(*l_iter) : section->getConstrainedFiberDimension(*l_iter, space);
368: }
369: }
370: if (debug) {std::cout << "["<<order->commRank()<<"] has offset " << val.prefix << " and size " << val.index << std::endl;}
371: localSize += val.index;
372: order->updatePoint(*l_iter, &val);
373: }
374: if (debug) {std::cout << "["<<order->commRank()<<"] local size" << localSize << std::endl;}
375: order->setLocalSize(localSize);
376: }
377: // Calculate process offsets
378: template<typename Numbering>
379: void calculateOffsets(const Obj<Numbering>& numbering) {
380: int localSize = numbering->getLocalSize();
381: int *offsets = new int[numbering->commSize()+1];
383: offsets[0] = 0;
384: MPI_Allgather(&localSize, 1, MPI_INT, &(offsets[1]), 1, MPI_INT, numbering->comm());
385: for(int p = 2; p <= numbering->commSize(); p++) {
386: offsets[p] += offsets[p-1];
387: }
388: numbering->setGlobalOffsets(offsets);
389: delete [] offsets;
390: }
391: // Update local offsets based upon process offsets
392: template<typename Numbering, typename Iterator>
393: void updateOrder(const Obj<Numbering>& numbering, const Iterator& pointsBegin, const Iterator& pointsEnd) {
394: const typename Numbering::value_type val = numbering->getGlobalOffset(numbering->commRank());
396: for(Iterator l_iter = pointsBegin; l_iter != pointsEnd; ++l_iter) {
397: if (numbering->isLocal(*l_iter)) {
398: numbering->updateAddPoint(*l_iter, &val);
399: }
400: }
401: }
402: // Communicate numbers in the overlap
403: void completeNumbering(const Obj<numbering_type>& numbering, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, bool allowDuplicates = false) {
404: #if 1
405: typedef ALE::UniformSection<ALE::Pair<int, typename send_overlap_type::source_type>, typename numbering_type::value_type> OverlapSection;
406: typedef typename OverlapSection::point_type overlap_point_type;
407: Obj<OverlapSection> overlapSection = new OverlapSection(numbering->comm(), numbering->debug());
408: const int debug = sendOverlap->debug();
410: if (debug) {std::cout << "["<<numbering->commRank()<<"] Completing numbering" << std::endl;}
411: ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, numbering, overlapSection);
412: if (debug) {overlapSection->view("Overlap Section");}
413: const typename recv_overlap_type::capSequence::iterator rBegin = recvOverlap->capBegin();
414: const typename recv_overlap_type::capSequence::iterator rEnd = recvOverlap->capEnd();
416: for(typename recv_overlap_type::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
417: const int rank = *r_iter;
418: const typename recv_overlap_type::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
419: const typename recv_overlap_type::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
421: for(typename recv_overlap_type::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
422: const point_type& localPoint = *p_iter;
423: const point_type& remotePoint = p_iter.color();
424: const overlap_point_type oPoint = overlap_point_type(rank, remotePoint);
425: const int fDim = overlapSection->getFiberDimension(oPoint);
426: const typename numbering_type::value_type *values = overlapSection->restrictPoint(oPoint);
428: for(int i = 0; i < fDim; ++i) {
429: if (debug) {std::cout << "["<<numbering->commRank()<<"] local point " << localPoint << " remote point " << remotePoint << " number " << values[i] << std::endl;}
430: if (values[i] >= 0) {
431: if (numbering->isLocal(localPoint) && !allowDuplicates) {
432: ostringstream msg;
433: msg << "["<<numbering->commRank()<<"]Multiple indices for local point " << localPoint << " remote point " << remotePoint << " from " << rank << " with index " << values[i];
434: throw ALE::Exception(msg.str().c_str());
435: }
436: if (!numbering->hasPoint(localPoint)) {
437: ostringstream msg;
438: msg << "["<<numbering->commRank()<<"]Unexpected local point " << localPoint << " remote point " << remotePoint << " from " << rank << " with index " << values[i];
439: throw ALE::Exception(msg.str().c_str());
440: }
441: int val = -(values[i]+1);
442: numbering->updatePoint(localPoint, &val);
443: }
444: }
445: }
446: }
447: #else
448: typedef Field<send_overlap_type, int, Section<point_type, value_type, value_alloc_type> > send_section_type;
449: typedef Field<recv_overlap_type, int, Section<point_type, value_type, value_alloc_type> > recv_section_type;
450: typedef typename ALE::DiscreteSieve<point_type, alloc_type> dsieve_type;
451: typedef typename ALE::Topology<int, dsieve_type, alloc_type> dtopology_type;
452: typedef typename ALE::New::SectionCompletion<dtopology_type, int, alloc_type> completion;
453: const Obj<send_section_type> sendSection = new send_section_type(numbering->comm(), this->debug());
454: const Obj<recv_section_type> recvSection = new recv_section_type(numbering->comm(), sendSection->getTag(), this->debug());
455: const int debug = sendOverlap->debug();
457: if (debug) {std::cout << "["<<numbering->commRank()<<"] Completing numbering" << std::endl;}
458: completion::completeSection(sendOverlap, recvOverlap, numbering->getAtlas(), numbering, sendSection, recvSection);
459: const Obj<typename recv_overlap_type::baseSequence> rPoints = recvOverlap->base();
461: for(typename recv_overlap_type::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rPoints->end(); ++p_iter) {
462: const Obj<typename recv_overlap_type::coneSequence>& ranks = recvOverlap->cone(*p_iter);
463: const typename recv_overlap_type::target_type& localPoint = *p_iter;
465: for(typename recv_overlap_type::coneSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
466: const typename recv_overlap_type::target_type& remotePoint = r_iter.color();
467: const int rank = *r_iter;
468: const Obj<typename recv_section_type::section_type>& section = recvSection->getSection(rank);
469: const typename recv_section_type::value_type *values = section->restrictPoint(remotePoint);
471: if (section->getFiberDimension(remotePoint) == 0) continue;
472: if (debug) {std::cout << "["<<numbering->commRank()<<"] local point " << localPoint << " remote point " << remotePoint << " number " << values[0] << std::endl;}
473: if (values[0] >= 0) {
474: if (debug) {std::cout << "["<<numbering->commRank()<<"] local point " << localPoint << " dim " << numbering->getAtlas()->getFiberDimension(localPoint) << std::endl;}
475: if (numbering->isLocal(localPoint) && !allowDuplicates) {
476: ostringstream msg;
477: msg << "["<<numbering->commRank()<<"]Multiple indices for local point " << localPoint << " remote point " << remotePoint << " from " << rank << " with index " << values[0];
478: throw ALE::Exception(msg.str().c_str());
479: }
480: if (numbering->getAtlas()->getFiberDimension(localPoint) == 0) {
481: ostringstream msg;
482: msg << "["<<numbering->commRank()<<"]Unexpected local point " << localPoint << " remote point " << remotePoint << " from " << rank << " with index " << values[0];
483: throw ALE::Exception(msg.str().c_str());
484: }
485: int val = -(values[0]+1);
486: numbering->updatePoint(localPoint, &val);
487: }
488: }
489: }
490: #endif
491: }
492: // Communicate (size,offset)s in the overlap
493: void completeOrder(const Obj<order_type>& order, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, bool allowDuplicates = false) {
494: #if 1
495: typedef ALE::UniformSection<ALE::Pair<int, typename send_overlap_type::source_type>, typename order_type::value_type> OverlapSection;
496: typedef typename OverlapSection::point_type overlap_point_type;
497: Obj<OverlapSection> overlapSection = new OverlapSection(order->comm(), order->debug());
498: const int debug = sendOverlap->debug();
500: if (debug) {std::cout << "["<<order->commRank()<<"] Completing ordering" << std::endl;}
501: ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, order, overlapSection);
502: if (debug) {overlapSection->view("Overlap Section");}
503: const typename recv_overlap_type::capSequence::iterator rBegin = recvOverlap->capBegin();
504: const typename recv_overlap_type::capSequence::iterator rEnd = recvOverlap->capEnd();
506: for(typename recv_overlap_type::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
507: const int rank = *r_iter;
508: const typename recv_overlap_type::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
509: const typename recv_overlap_type::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
511: for(typename recv_overlap_type::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
512: const point_type& localPoint = *p_iter;
513: const point_type& remotePoint = p_iter.color();
514: const overlap_point_type oPoint = overlap_point_type(rank, remotePoint);
515: const int fDim = overlapSection->getFiberDimension(oPoint);
516: const typename order_type::value_type *values = overlapSection->restrictPoint(oPoint);
518: for(int i = 0; i < fDim; ++i) {
519: if (debug) {std::cout << "["<<order->commRank()<<"] local point " << localPoint << " remote point " << remotePoint<<"("<<rank<<")" << " offset " << values[i].prefix << " and size " << values[i].index << std::endl;}
520: if (values[i].index == 0) continue;
521: if (values[i].prefix >= 0) {
522: if (order->isLocal(localPoint)) {
523: if (!allowDuplicates) {
524: ostringstream msg;
525: msg << "["<<order->commRank()<<"]Multiple indices for local point " << localPoint << " remote point " << remotePoint << " from " << rank << " with index " << values[i];
526: throw ALE::Exception(msg.str().c_str());
527: }
528: continue;
529: }
530: const oValue_type val(-(values[i].prefix+1), values[i].index);
531: order->updatePoint(localPoint, &val);
532: } else {
533: if (order->isLocal(localPoint)) continue;
534: order->updatePoint(localPoint, &values[i]);
535: }
536: }
537: }
538: }
539: #else
540: typedef Field<send_overlap_type, int, Section<point_type, oValue_type, oValue_alloc_type> > send_section_type;
541: typedef Field<recv_overlap_type, int, Section<point_type, oValue_type, oValue_alloc_type> > recv_section_type;
542: typedef ConstantSection<point_type, int, alloc_type> constant_sizer;
543: typedef typename ALE::DiscreteSieve<point_type, alloc_type> dsieve_type;
544: typedef typename ALE::Topology<int, dsieve_type, alloc_type> dtopology_type;
545: typedef typename ALE::New::SectionCompletion<dtopology_type, int, alloc_type> completion;
546: const Obj<send_section_type> sendSection = new send_section_type(order->comm(), this->debug());
547: const Obj<recv_section_type> recvSection = new recv_section_type(order->comm(), sendSection->getTag(), this->debug());
548: const int debug = sendOverlap->debug();
550: if (debug) {std::cout << "["<<order->commRank()<<"] Completing ordering" << std::endl;}
551: completion::completeSection(sendOverlap, recvOverlap, order->getAtlas(), order, sendSection, recvSection);
552: Obj<typename recv_overlap_type::baseSequence> recvPoints = recvOverlap->base();
554: for(typename recv_overlap_type::baseSequence::iterator p_iter = recvPoints->begin(); p_iter != recvPoints->end(); ++p_iter) {
555: if (!order->hasPoint(*p_iter)) {
556: order->setFiberDimension(*p_iter, 1);
557: order->updatePoint(*p_iter, &this->_unknownOrder);
558: }
559: }
560: for(typename recv_overlap_type::baseSequence::iterator p_iter = recvPoints->begin(); p_iter != recvPoints->end(); ++p_iter) {
561: const Obj<typename recv_overlap_type::coneSequence>& ranks = recvOverlap->cone(*p_iter);
562: const typename recv_overlap_type::target_type& localPoint = *p_iter;
564: for(typename recv_overlap_type::coneSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
565: const typename recv_overlap_type::target_type& remotePoint = r_iter.color();
566: const int rank = *r_iter;
567: const Obj<typename recv_section_type::section_type>& section = recvSection->getSection(rank);
568: const typename recv_section_type::value_type *values = section->restrictPoint(remotePoint);
570: if (section->getFiberDimension(remotePoint) == 0) continue;
571: if (debug) {std::cout << "["<<order->commRank()<<"] local point " << localPoint << " remote point " << remotePoint<<"("<<rank<<")" << " offset " << values[0].prefix << " and size " << values[0].index << std::endl;}
572: if (values[0].index == 0) continue;
573: if (values[0].prefix >= 0) {
574: if (order->isLocal(localPoint)) {
575: if (!allowDuplicates) {
576: ostringstream msg;
577: msg << "["<<order->commRank()<<"]Multiple indices for local point " << localPoint << " remote point " << remotePoint << " from " << rank << " with index " << values[0];
578: throw ALE::Exception(msg.str().c_str());
579: }
580: continue;
581: }
582: const oValue_type val(-(values[0].prefix+1), values[0].index);
583: order->updatePoint(localPoint, &val);
584: } else {
585: if (order->isLocal(localPoint)) continue;
586: order->updatePoint(localPoint, values);
587: }
588: }
589: }
590: #endif
591: }
592: // Construct a full global numbering
593: template<typename Iterator>
594: void constructNumbering(const Obj<numbering_type>& numbering, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Iterator& pointsBegin, const Iterator& pointsEnd) {
595: this->constructLocalNumbering(numbering, sendOverlap, pointsBegin, pointsEnd);
596: this->calculateOffsets(numbering);
597: this->updateOrder(numbering, pointsBegin, pointsEnd);
598: this->completeNumbering(numbering, sendOverlap, recvOverlap);
599: }
600: // Construct a full global order
601: template<typename Iterator, typename Section>
602: void constructOrder(const Obj<order_type>& order, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Iterator& pointsBegin, const Iterator& pointsEnd, const Obj<Section>& section, const int space = -1, const Obj<label_type>& label = PETSC_NULL) {
603: this->constructLocalOrder(order, sendOverlap, pointsBegin, pointsEnd, section, space, false, label);
604: this->calculateOffsets(order);
605: this->updateOrder(order, pointsBegin, pointsEnd);
606: this->completeOrder(order, sendOverlap, recvOverlap);
607: }
608: template<typename Iterator, typename Section>
609: void constructOrderBC(const Obj<order_type>& order, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Iterator& pointsBegin, const Iterator& pointsEnd, const Obj<Section>& section, const int space = -1, const Obj<label_type>& label = PETSC_NULL) {
610: this->constructLocalOrder(order, sendOverlap, pointsBegin, pointsEnd, section, space, true, label);
611: this->calculateOffsets(order);
612: this->updateOrder(order, pointsBegin, pointsEnd);
613: this->completeOrder(order, sendOverlap, recvOverlap);
614: }
615: public:
616: // Construct the inverse map from numbers to points
617: // If we really need this, then we should consider using a label
618: void constructInverseOrder(const Obj<numbering_type>& numbering) {
619: const typename numbering_type::chart_type& chart = numbering->getChart();
621: for(typename numbering_type::chart_type::iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
622: numbering->setPoint(numbering->getIndex(*p_iter), *p_iter);
623: }
624: }
625: public: // Real interface
626: template<typename ABundle_>
627: const Obj<numbering_type>& getLocalNumbering(const Obj<ABundle_>& bundle, const int depth) {
628: if ((this->_localNumberings.find(bundle.ptr()) == this->_localNumberings.end()) ||
629: (this->_localNumberings[bundle.ptr()].find("depth") == this->_localNumberings[bundle.ptr()].end()) ||
630: (this->_localNumberings[bundle.ptr()]["depth"].find(depth) == this->_localNumberings[bundle.ptr()]["depth"].end())) {
631: Obj<numbering_type> numbering = new numbering_type(bundle->comm(), bundle->debug());
632: Obj<send_overlap_type> sendOverlap = new send_overlap_type(bundle->comm(), bundle->debug());
634: this->constructLocalNumbering(numbering, sendOverlap, bundle->depthStratum(depth)->begin(), bundle->depthStratum(depth)->end());
635: if (this->_debug) {std::cout << "Creating new local numbering: ptr " << bundle.ptr() << " depth " << depth << std::endl;}
636: this->_localNumberings[bundle.ptr()]["depth"][depth] = numbering;
637: } else {
638: if (this->_debug) {std::cout << "Using old local numbering: ptr " << bundle.ptr() << " depth " << depth << std::endl;}
639: }
640: return this->_localNumberings[bundle.ptr()]["depth"][depth];
641: }
642: template<typename ABundle_>
643: const Obj<numbering_type>& getNumbering(const Obj<ABundle_>& bundle, const int depth) {
644: if ((this->_numberings.find(bundle.ptr()) == this->_numberings.end()) ||
645: (this->_numberings[bundle.ptr()].find("depth") == this->_numberings[bundle.ptr()].end()) ||
646: (this->_numberings[bundle.ptr()]["depth"].find(depth) == this->_numberings[bundle.ptr()]["depth"].end())) {
647: bundle->constructOverlap();
648: Obj<numbering_type> numbering = new numbering_type(bundle->comm(), bundle->debug());
649: Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
650: Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();
652: if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Creating new numbering: fixed depth value " << depth << std::endl;}
653: if (depth == -1) {
654: this->constructNumbering(numbering, sendOverlap, recvOverlap, bundle->getSieve()->getChart().begin(), bundle->getSieve()->getChart().end());
655: } else {
656: this->constructNumbering(numbering, sendOverlap, recvOverlap, bundle->depthStratum(depth)->begin(), bundle->depthStratum(depth)->end());
657: }
658: this->_numberings[bundle.ptr()]["depth"][depth] = numbering;
659: } else {
660: if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Using old numbering: fixed depth value " << depth << std::endl;}
661: }
662: return this->_numberings[bundle.ptr()]["depth"][depth];
663: }
664: template<typename ABundle_>
665: const Obj<numbering_type>& getNumbering(const Obj<ABundle_>& bundle, const std::string& labelname, const int value) {
666: if ((this->_numberings.find(bundle.ptr()) == this->_numberings.end()) ||
667: (this->_numberings[bundle.ptr()].find(labelname) == this->_numberings[bundle.ptr()].end()) ||
668: (this->_numberings[bundle.ptr()][labelname].find(value) == this->_numberings[bundle.ptr()][labelname].end())) {
669: bundle->constructOverlap();
670: Obj<numbering_type> numbering = new numbering_type(bundle->comm(), bundle->debug());
671: Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
672: Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();
674: if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Creating new numbering: " << labelname << " value " << value << std::endl;}
675: this->constructNumbering(numbering, sendOverlap, recvOverlap, bundle->getLabelStratum(labelname, value)->begin(), bundle->getLabelStratum(labelname, value)->end());
676: this->_numberings[bundle.ptr()][labelname][value] = numbering;
677: } else {
678: if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Using old numbering: " << labelname << " value " << value << std::endl;}
679: }
680: return this->_numberings[bundle.ptr()][labelname][value];
681: }
682: template<typename ABundle_, typename Section_>
683: const Obj<order_type>& getLocalOrder(const Obj<ABundle_>& bundle, const std::string& name, const Obj<Section_>& section, const int space = -1, const Obj<label_type>& label = PETSC_NULL) {
684: if ((this->_localOrders.find(bundle.ptr()) == this->_localOrders.end()) ||
685: (this->_localOrders[bundle.ptr()].find(name) == this->_localOrders[bundle.ptr()].end())) {
686: Obj<order_type> order = new order_type(bundle->comm(), bundle->debug());
687: Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
688: Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();
690: if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Creating new local order: " << name << std::endl;}
691: this->constructLocalOrder(order, sendOverlap, section->getChart().begin(), section->getChart().end(), section, space, false, label);
692: this->_localOrders[bundle.ptr()][name] = order;
693: } else {
694: if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Using old local order: " << name << std::endl;}
695: }
696: return this->_localOrders[bundle.ptr()][name];
697: }
698: template<typename ABundle_, typename Section_>
699: const Obj<order_type>& getGlobalOrder(const Obj<ABundle_>& bundle, const std::string& name, const Obj<Section_>& section, const int space = -1, const Obj<label_type>& label = PETSC_NULL) {
700: if ((this->_orders.find(bundle.ptr()) == this->_orders.end()) ||
701: (this->_orders[bundle.ptr()].find(name) == this->_orders[bundle.ptr()].end())) {
702: bundle->constructOverlap();
703: Obj<order_type> order = new order_type(bundle->comm(), bundle->debug());
704: Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
705: Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();
707: if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Creating new global order: " << name << std::endl;}
708: this->constructOrder(order, sendOverlap, recvOverlap, section->getChart().begin(), section->getChart().end(), section, space, label);
709: this->_orders[bundle.ptr()][name] = order;
710: } else {
711: if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Using old global order: " << name << std::endl;}
712: }
713: return this->_orders[bundle.ptr()][name];
714: }
715: template<typename ABundle_, typename Iterator_, typename Section_>
716: const Obj<order_type>& getGlobalOrder(const Obj<ABundle_>& bundle, const std::string& name, const Iterator_& pointsBegin, const Iterator_& pointsEnd, const Obj<Section_>& section, const int space = -1, const Obj<label_type>& label = PETSC_NULL) {
717: if ((this->_orders.find(bundle.ptr()) == this->_orders.end()) ||
718: (this->_orders[bundle.ptr()].find(name) == this->_orders[bundle.ptr()].end())) {
719: bundle->constructOverlap();
720: Obj<order_type> order = new order_type(bundle->comm(), bundle->debug());
721: Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
722: Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();
724: if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Creating new global order: " << name << std::endl;}
725: this->constructOrder(order, sendOverlap, recvOverlap, pointsBegin, pointsEnd, section, space, label);
726: this->_orders[bundle.ptr()][name] = order;
727: } else {
728: if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Using old global order: " << name << std::endl;}
729: }
730: return this->_orders[bundle.ptr()][name];
731: }
732: template<typename ABundle_, typename Section_>
733: const Obj<order_type>& getGlobalOrderWithBC(const Obj<ABundle_>& bundle, const std::string& name, const Obj<Section_>& section, const int space = -1, const Obj<label_type>& label = PETSC_NULL) {
734: if ((this->_orders.find(bundle.ptr()) == this->_orders.end()) ||
735: (this->_orders[bundle.ptr()].find(name) == this->_orders[bundle.ptr()].end())) {
736: bundle->constructOverlap();
737: Obj<order_type> order = new order_type(bundle->comm(), bundle->debug());
738: Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
739: Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();
741: if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Creating new global order: " << name << std::endl;}
742: this->constructOrderBC(order, sendOverlap, recvOverlap, section->getChart().begin(), section->getChart().end(), section, space, label);
743: this->_orders[bundle.ptr()][name] = order;
744: } else {
745: if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Using old global order: " << name << std::endl;}
746: }
747: return this->_orders[bundle.ptr()][name];
748: }
749: template<typename ABundle_, typename Iterator_, typename Section_>
750: const Obj<order_type>& getGlobalOrderWithBC(const Obj<ABundle_>& bundle, const std::string& name, const Iterator_& pointsBegin, const Iterator_& pointsEnd, const Obj<Section_>& section, const int space = -1, const Obj<label_type>& label = PETSC_NULL) {
751: if ((this->_ordersBC.find(bundle.ptr()) == this->_ordersBC.end()) ||
752: (this->_ordersBC[bundle.ptr()].find(name) == this->_ordersBC[bundle.ptr()].end())) {
753: bundle->constructOverlap();
754: Obj<order_type> order = new order_type(bundle->comm(), bundle->debug());
755: Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
756: Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();
758: if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Creating new global order: " << name << std::endl;}
759: this->constructOrderBC(order, sendOverlap, recvOverlap, pointsBegin, pointsEnd, section, space, label);
760: this->_orders[bundle.ptr()][name] = order;
761: } else {
762: if (this->_debug) {std::cout << "["<<bundle->commRank()<<"]Using old global order: " << name << std::endl;}
763: }
764: return this->_orders[bundle.ptr()][name];
765: }
766: template<typename ABundle_>
767: void setGlobalOrder(const Obj<ABundle_>& bundle, const std::string& name, const Obj<order_type>& order) {
768: this->_orders[bundle.ptr()][name] = order;
769: }
770: };
771: }
772: #endif