Actual source code: Topology.hh
petsc-3.3-p7 2013-05-11
1: #ifndef included_ALE_CoSieve_hh
2: #define included_ALE_CoSieve_hh
4: #ifndef included_ALE_Sieve_hh
5: #include <sieve/Sieve.hh>
6: #endif
8: namespace ALE {
9: // A Topology is a collection of Sieves
10: // Each Sieve has a label, which we call a \emph{patch}
11: // The collection itself we call a \emph{sheaf}
12: // The main operation we provide in Topology is the creation of a \emph{label}
13: // A label is a bidirectional mapping of Sieve points to integers, implemented with a Sifter
14: template<typename Patch_, typename Sieve_, typename Alloc_ = malloc_allocator<typename Sieve_::point_type> >
15: class Topology : public ALE::ParallelObject {
16: public:
17: typedef Patch_ patch_type;
18: typedef Sieve_ sieve_type;
19: typedef Alloc_ alloc_type;
20: typedef typename alloc_type::template rebind<int>::other int_alloc_type;
21: typedef typename sieve_type::point_type point_type;
22: typedef typename alloc_type::template rebind<std::pair<const patch_type, Obj<sieve_type> > >::other pair1_alloc_type;
23: typedef typename std::map<patch_type, Obj<sieve_type>, std::less<patch_type>, pair1_alloc_type> sheaf_type;
24: typedef typename ALE::Sifter<int, point_type, int> patch_label_type;
25: typedef typename alloc_type::template rebind<patch_label_type>::other patch_label_alloc_type;
26: typedef typename patch_label_alloc_type::pointer patch_label_ptr;
27: typedef typename alloc_type::template rebind<std::pair<const patch_type, Obj<patch_label_type> > >::other pair2_alloc_type;
28: typedef typename std::map<patch_type, Obj<patch_label_type>, std::less<patch_type>, pair2_alloc_type> label_type;
29: typedef typename alloc_type::template rebind<std::pair<const patch_type, int> >::other pair3_alloc_type;
30: typedef typename std::map<patch_type, int, std::less<patch_type>, pair3_alloc_type> max_label_type;
31: typedef typename alloc_type::template rebind<std::pair<const std::string, label_type> >::other pair4_alloc_type;
32: typedef typename std::map<const std::string, label_type, std::less<const std::string>, pair4_alloc_type> labels_type;
33: typedef typename patch_label_type::supportSequence label_sequence;
34: typedef typename std::set<point_type, std::less<point_type>, alloc_type> point_set_type;
35: typedef typename alloc_type::template rebind<point_set_type>::other point_set_alloc_type;
36: typedef typename point_set_alloc_type::pointer point_set_ptr;
37: typedef typename ALE::Sifter<int,point_type,point_type> send_overlap_type;
38: typedef typename alloc_type::template rebind<send_overlap_type>::other send_overlap_alloc_type;
39: typedef typename send_overlap_alloc_type::pointer send_overlap_ptr;
40: typedef typename ALE::Sifter<point_type,int,point_type> recv_overlap_type;
41: typedef typename alloc_type::template rebind<recv_overlap_type>::other recv_overlap_alloc_type;
42: typedef typename recv_overlap_alloc_type::pointer recv_overlap_ptr;
43: protected:
44: sheaf_type _sheaf;
45: labels_type _labels;
46: int _maxHeight;
47: max_label_type _maxHeights;
48: int _maxDepth;
49: max_label_type _maxDepths;
50: bool _calculatedOverlap;
51: Obj<send_overlap_type> _sendOverlap;
52: Obj<recv_overlap_type> _recvOverlap;
53: Obj<send_overlap_type> _distSendOverlap;
54: Obj<recv_overlap_type> _distRecvOverlap;
55: // Work space
56: Obj<point_set_type> _modifiedPoints;
57: alloc_type _allocator;
58: int_alloc_type _int_allocator;
59: public:
60: Topology(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug), _maxHeight(-1), _maxDepth(-1), _calculatedOverlap(false) {
61: send_overlap_ptr pSendOverlap = send_overlap_alloc_type(this->_allocator).allocate(1);
62: send_overlap_alloc_type(this->_allocator).construct(pSendOverlap, send_overlap_type(this->comm(), this->debug()));
63: this->_sendOverlap = Obj<send_overlap_type>(pSendOverlap, sizeof(send_overlap_type));
64: ///this->_sendOverlap = new send_overlap_type(this->comm(), this->debug());
65: recv_overlap_ptr pRecvOverlap = recv_overlap_alloc_type(this->_allocator).allocate(1);
66: recv_overlap_alloc_type(this->_allocator).construct(pRecvOverlap, recv_overlap_type(this->comm(), this->debug()));
67: this->_recvOverlap = Obj<recv_overlap_type>(pRecvOverlap, sizeof(recv_overlap_type));
68: ///this->_recvOverlap = new recv_overlap_type(this->comm(), this->debug());
69: point_set_ptr pModPoints = point_set_alloc_type(this->_allocator).allocate(1);
70: point_set_alloc_type(this->_allocator).construct(pModPoints, point_set_type());
71: this->_modifiedPoints = Obj<point_set_type>(pModPoints, sizeof(point_set_type));
72: ///this->_modifiedPoints = new point_set_type();
73: };
74: virtual ~Topology() {};
75: public: // Verifiers
76: void checkPatch(const patch_type& patch) {
77: if (this->_sheaf.find(patch) == this->_sheaf.end()) {
78: ostringstream msg;
79: msg << "Invalid topology patch: " << patch << std::endl;
80: throw ALE::Exception(msg.str().c_str());
81: }
82: };
83: void checkLabel(const std::string& name, const patch_type& patch) {
84: this->checkPatch(patch);
85: if ((this->_labels.find(name) == this->_labels.end()) || (this->_labels[name].find(patch) == this->_labels[name].end())) {
86: ostringstream msg;
87: msg << "Invalid label name: " << name << " for patch " << patch << std::endl;
88: throw ALE::Exception(msg.str().c_str());
89: }
90: };
91: bool hasPatch(const patch_type& patch) {
92: if (this->_sheaf.find(patch) != this->_sheaf.end()) {
93: return true;
94: }
95: return false;
96: };
97: bool hasLabel(const std::string& name, const patch_type& patch) {
98: if ((this->_labels.find(name) != this->_labels.end()) && (this->_labels[name].find(patch) != this->_labels[name].end())) {
99: return true;
100: }
101: return false;
102: };
103: public: // Accessors
104: const Obj<sieve_type>& getPatch(const patch_type& patch) {
105: this->checkPatch(patch);
106: return this->_sheaf[patch];
107: };
108: void setPatch(const patch_type& patch, const Obj<sieve_type>& sieve) {
109: this->_sheaf[patch] = sieve;
110: };
111: int getValue (const Obj<patch_label_type>& label, const point_type& point, const int defValue = 0) {
112: const Obj<typename patch_label_type::coneSequence>& cone = label->cone(point);
114: if (cone->size() == 0) return defValue;
115: return *cone->begin();
116: };
117: template<typename InputPoints>
118: int getMaxValue (const Obj<patch_label_type>& label, const Obj<InputPoints>& points, const int defValue = 0) {
119: int maxValue = defValue;
121: for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
122: maxValue = std::max(maxValue, this->getValue(label, *p_iter, defValue));
123: }
124: return maxValue;
125: }
126: void setValue(const Obj<patch_label_type>& label, const point_type& point, const int value) {
127: label->setCone(value, point);
128: };
129: const Obj<patch_label_type>& createLabel(const patch_type& patch, const std::string& name) {
130: this->checkPatch(patch);
131: ///patch_label_ptr pLabel = patch_label_alloc_type(this->_allocator).allocate(1);
132: ///patch_label_alloc_type(this->_allocator).construct(pLabel, patch_label_type(this->comm(), this->debug()));
133: ///this->_labels[name][patch] = Obj<patch_label_type>(pLabel, sizeof(patch_label_type));
134: this->_labels[name][patch] = new patch_label_type(this->comm(), this->debug());
135: return this->_labels[name][patch];
136: };
137: const Obj<patch_label_type>& getLabel(const patch_type& patch, const std::string& name) {
138: this->checkLabel(name, patch);
139: return this->_labels[name][patch];
140: };
141: const Obj<label_sequence>& getLabelStratum(const patch_type& patch, const std::string& name, int value) {
142: this->checkLabel(name, patch);
143: return this->_labels[name][patch]->support(value);
144: };
145: const sheaf_type& getPatches() {
146: return this->_sheaf;
147: };
148: const labels_type& getLabels() {
149: return this->_labels;
150: };
151: void clear() {
152: this->_sheaf.clear();
153: this->_labels.clear();
154: this->_maxHeight = -1;
155: this->_maxHeights.clear();
156: this->_maxDepth = -1;
157: this->_maxDepths.clear();
158: };
159: const Obj<send_overlap_type>& getSendOverlap() const {return this->_sendOverlap;};
160: void setSendOverlap(const Obj<send_overlap_type>& overlap) {this->_sendOverlap = overlap;};
161: const Obj<recv_overlap_type>& getRecvOverlap() const {return this->_recvOverlap;};
162: void setRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_recvOverlap = overlap;};
163: const Obj<send_overlap_type>& getDistSendOverlap() const {return this->_distSendOverlap;};
164: void setDistSendOverlap(const Obj<send_overlap_type>& overlap) {this->_distSendOverlap = overlap;};
165: const Obj<recv_overlap_type>& getDistRecvOverlap() const {return this->_distRecvOverlap;};
166: void setDistRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_distRecvOverlap = overlap;};
167: public: // Stratification
168: template<class InputPoints>
169: void computeHeight(const Obj<patch_label_type>& height, const Obj<sieve_type>& sieve, const Obj<InputPoints>& points, int& maxHeight) {
170: this->_modifiedPoints->clear();
172: for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
173: // Compute the max height of the points in the support of p, and add 1
174: int h0 = this->getValue(height, *p_iter, -1);
175: int h1 = this->getMaxValue(height, sieve->support(*p_iter), -1) + 1;
177: if(h1 != h0) {
178: this->setValue(height, *p_iter, h1);
179: if (h1 > maxHeight) maxHeight = h1;
180: this->_modifiedPoints->insert(*p_iter);
181: }
182: }
183: // FIX: We would like to avoid the copy here with cone()
184: if(this->_modifiedPoints->size() > 0) {
185: this->computeHeight(height, sieve, sieve->cone(this->_modifiedPoints), maxHeight);
186: }
187: }
188: void computeHeights() {
189: const std::string name("height");
191: this->_maxHeight = -1;
192: for(typename sheaf_type::iterator s_iter = this->_sheaf.begin(); s_iter != this->_sheaf.end(); ++s_iter) {
193: const Obj<patch_label_type>& label = this->createLabel(s_iter->first, name);
195: this->_maxHeights[s_iter->first] = -1;
196: this->computeHeight(label, s_iter->second, s_iter->second->leaves(), this->_maxHeights[s_iter->first]);
197: if (this->_maxHeights[s_iter->first] > this->_maxHeight) this->_maxHeight = this->_maxHeights[s_iter->first];
198: }
199: }
200: int height() const {return this->_maxHeight;};
201: int height(const patch_type& patch) {
202: this->checkPatch(patch);
203: return this->_maxHeights[patch];
204: }
205: int height(const patch_type& patch, const point_type& point) {
206: return this->getValue(this->_labels["height"][patch], point, -1);
207: }
208: const Obj<label_sequence>& heightStratum(const patch_type& patch, int height) {
209: return this->getLabelStratum(patch, "height", height);
210: }
211: template<class InputPoints>
212: void computeDepth(const Obj<patch_label_type>& depth, const Obj<sieve_type>& sieve, const Obj<InputPoints>& points, int& maxDepth) {
213: this->_modifiedPoints->clear();
215: for(typename InputPoints::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
216: // Compute the max depth of the points in the cone of p, and add 1
217: int d0 = this->getValue(depth, *p_iter, -1);
218: int d1 = this->getMaxValue(depth, sieve->cone(*p_iter), -1) + 1;
220: if(d1 != d0) {
221: this->setValue(depth, *p_iter, d1);
222: if (d1 > maxDepth) maxDepth = d1;
223: this->_modifiedPoints->insert(*p_iter);
224: }
225: }
226: // FIX: We would like to avoid the copy here with support()
227: if(this->_modifiedPoints->size() > 0) {
228: this->computeDepth(depth, sieve, sieve->support(this->_modifiedPoints), maxDepth);
229: }
230: }
231: void computeDepths() {
232: const std::string name("depth");
234: this->_maxDepth = -1;
235: for(typename sheaf_type::iterator s_iter = this->_sheaf.begin(); s_iter != this->_sheaf.end(); ++s_iter) {
236: const Obj<patch_label_type>& label = this->createLabel(s_iter->first, name);
238: this->_maxDepths[s_iter->first] = -1;
239: this->computeDepth(label, s_iter->second, s_iter->second->roots(), this->_maxDepths[s_iter->first]);
240: if (this->_maxDepths[s_iter->first] > this->_maxDepth) this->_maxDepth = this->_maxDepths[s_iter->first];
241: }
242: }
243: int depth() const {return this->_maxDepth;};
244: int depth(const patch_type& patch) {
245: this->checkPatch(patch);
246: return this->_maxDepths[patch];
247: }
248: int depth(const patch_type& patch, const point_type& point) {
249: return this->getValue(this->_labels["depth"][patch], point, -1);
250: }
251: const Obj<label_sequence>& depthStratum(const patch_type& patch, int depth) {
252: return this->getLabelStratum(patch, "depth", depth);
253: }
256: void stratify() {
257: ALE_LOG_EVENT_BEGIN;
258: this->computeHeights();
259: this->computeDepths();
260: ALE_LOG_EVENT_END;
261: }
262: public: // Viewers
263: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
264: if (comm == MPI_COMM_NULL) {
265: comm = this->comm();
266: }
267: if (name == "") {
268: PetscPrintf(comm, "viewing a Topology\n");
269: } else {
270: PetscPrintf(comm, "viewing Topology '%s'\n", name.c_str());
271: }
272: PetscPrintf(comm, " maximum height %d maximum depth %d\n", this->height(), this->depth());
273: for(typename sheaf_type::const_iterator s_iter = this->_sheaf.begin(); s_iter != this->_sheaf.end(); ++s_iter) {
274: ostringstream txt;
276: txt << "Patch " << s_iter->first;
277: s_iter->second->view(txt.str().c_str(), comm);
278: PetscPrintf(comm, " maximum height %d maximum depth %d\n", this->height(s_iter->first), this->depth(s_iter->first));
279: }
280: for(typename labels_type::const_iterator l_iter = this->_labels.begin(); l_iter != this->_labels.end(); ++l_iter) {
281: PetscPrintf(comm, " label %s constructed\n", l_iter->first.c_str());
282: }
283: }
284: public:
285: void constructOverlap(const patch_type& patch) {
286: if (this->_calculatedOverlap) return;
287: if (this->hasPatch(patch)) {
288: this->constructOverlap(this->getPatch(patch)->base(), this->_sendOverlap, this->_recvOverlap);
289: this->constructOverlap(this->getPatch(patch)->cap(), this->_sendOverlap, this->_recvOverlap);
290: }
291: if (this->debug()) {
292: this->_sendOverlap->view("Send overlap");
293: this->_recvOverlap->view("Receive overlap");
294: }
295: this->_calculatedOverlap = true;
296: }
297: template<typename Sequence>
298: void constructOverlap(const Obj<Sequence>& points, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
299: point_type *sendBuf = this->_allocator.allocate(points->size());
300: for(unsigned int i = 0; i < points->size(); ++i) {this->_allocator.construct(sendBuf+i, point_type());}
301: ///point_type *sendBuf = new point_type[points->size()];
302: int size = 0;
303: for(typename Sequence::iterator l_iter = points->begin(); l_iter != points->end(); ++l_iter) {
304: sendBuf[size++] = *l_iter;
305: }
306: int *sizes = this->_int_allocator.allocate(this->commSize()+1); // The number of points coming from each process
307: int *offsets = this->_int_allocator.allocate(this->commSize()+1); // Prefix sums for sizes
308: int *oldOffs = this->_int_allocator.allocate(this->commSize()+1); // Temporary storage
309: for(int i = 0; i < this->commSize()+1; ++i) {
310: this->_int_allocator.construct(sizes+i, 0);
311: this->_int_allocator.construct(offsets+i, 0);
312: this->_int_allocator.construct(oldOffs+i, 0);
313: }
314: ///int *sizes = new int[this->commSize()];
315: ///int *offsets = new int[this->commSize()+1];
316: ///int *oldOffs = new int[this->commSize()+1];
317: point_type *remotePoints = NULL; // The points from each process
318: int *remoteRanks = NULL; // The rank and number of overlap points of each process that overlaps another
320: // Change to Allgather() for the correct binning algorithm
321: MPI_Gather(&size, 1, MPI_INT, sizes, 1, MPI_INT, 0, this->comm());
322: if (this->commRank() == 0) {
323: offsets[0] = 0;
324: for(int p = 1; p <= this->commSize(); p++) {
325: offsets[p] = offsets[p-1] + sizes[p-1];
326: }
327: remotePoints = this->_allocator.allocate(offsets[this->commSize()]);
328: for(int i = 0; i < offsets[this->commSize()]; ++i) {this->_allocator.construct(remotePoints+i, point_type());}
329: ///remotePoints = new point_type[offsets[this->commSize()]];
330: }
331: MPI_Gatherv(sendBuf, size, MPI_INT, remotePoints, sizes, offsets, MPI_INT, 0, this->comm());
332: std::map<int, std::map<int, std::set<point_type> > > overlapInfo; // Maps (p,q) to their set of overlap points
334: if (this->commRank() == 0) {
335: for(int p = 0; p < this->commSize(); p++) {
336: std::sort(&remotePoints[offsets[p]], &remotePoints[offsets[p+1]]);
337: }
338: for(int p = 0; p <= this->commSize(); p++) {
339: oldOffs[p] = offsets[p];
340: }
341: for(int p = 0; p < this->commSize(); p++) {
342: for(int q = p+1; q < this->commSize(); q++) {
343: std::set_intersection(&remotePoints[oldOffs[p]], &remotePoints[oldOffs[p+1]],
344: &remotePoints[oldOffs[q]], &remotePoints[oldOffs[q+1]],
345: std::insert_iterator<std::set<point_type> >(overlapInfo[p][q], overlapInfo[p][q].begin()));
346: overlapInfo[q][p] = overlapInfo[p][q];
347: }
348: sizes[p] = overlapInfo[p].size()*2;
349: offsets[p+1] = offsets[p] + sizes[p];
350: }
351: remoteRanks = this->_int_allocator.allocate(offsets[this->commSize()]);
352: for(int i = 0; i < offsets[this->commSize()]; ++i) {this->_int_allocator.construct(remoteRanks+i, 0);}
353: ///remoteRanks = new int[offsets[this->commSize()]];
354: int k = 0;
355: for(int p = 0; p < this->commSize(); p++) {
356: for(typename std::map<int, std::set<point_type> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
357: remoteRanks[k*2] = r_iter->first;
358: remoteRanks[k*2+1] = r_iter->second.size();
359: k++;
360: }
361: }
362: }
363: int numOverlaps; // The number of processes overlapping this process
364: MPI_Scatter(sizes, 1, MPI_INT, &numOverlaps, 1, MPI_INT, 0, this->comm());
365: int *overlapRanks = this->_int_allocator.allocate(numOverlaps);
366: for(int i = 0; i < numOverlaps; ++i) {this->_int_allocator.construct(overlapRanks+i, 0);}
367: ///int *overlapRanks = new int[numOverlaps]; // The rank and overlap size for each overlapping process
368: MPI_Scatterv(remoteRanks, sizes, offsets, MPI_INT, overlapRanks, numOverlaps, MPI_INT, 0, this->comm());
369: point_type *sendPoints = NULL; // The points to send to each process
370: if (this->commRank() == 0) {
371: for(int p = 0, k = 0; p < this->commSize(); p++) {
372: sizes[p] = 0;
373: for(int r = 0; r < (int) overlapInfo[p].size(); r++) {
374: sizes[p] += remoteRanks[k*2+1];
375: k++;
376: }
377: offsets[p+1] = offsets[p] + sizes[p];
378: }
379: sendPoints = this->_allocator.allocate(offsets[this->commSize()]);
380: for(int i = 0; i < offsets[this->commSize()]; ++i) {this->_allocator.construct(sendPoints+i, point_type());}
381: ///sendPoints = new point_type[offsets[this->commSize()]];
382: for(int p = 0, k = 0; p < this->commSize(); p++) {
383: for(typename std::map<int, std::set<point_type> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
384: int rank = r_iter->first;
385: for(typename std::set<point_type>::iterator p_iter = (overlapInfo[p][rank]).begin(); p_iter != (overlapInfo[p][rank]).end(); ++p_iter) {
386: sendPoints[k++] = *p_iter;
387: }
388: }
389: }
390: }
391: int numOverlapPoints = 0;
392: for(int r = 0; r < numOverlaps/2; r++) {
393: numOverlapPoints += overlapRanks[r*2+1];
394: }
395: point_type *overlapPoints = this->_allocator.allocate(numOverlapPoints);
396: for(int i = 0; i < numOverlapPoints; ++i) {this->_allocator.construct(overlapPoints+i, point_type());}
397: ///point_type *overlapPoints = new point_type[numOverlapPoints];
398: MPI_Scatterv(sendPoints, sizes, offsets, MPI_INT, overlapPoints, numOverlapPoints, MPI_INT, 0, this->comm());
400: for(int r = 0, k = 0; r < numOverlaps/2; r++) {
401: int rank = overlapRanks[r*2];
403: for(int p = 0; p < overlapRanks[r*2+1]; p++) {
404: point_type point = overlapPoints[k++];
406: sendOverlap->addArrow(point, rank, point);
407: recvOverlap->addArrow(rank, point, point);
408: }
409: }
411: delete [] overlapPoints;
412: delete [] overlapRanks;
413: delete [] sizes;
414: delete [] offsets;
415: delete [] oldOffs;
416: if (this->commRank() == 0) {
417: delete [] remoteRanks;
418: delete [] remotePoints;
419: delete [] sendPoints;
420: }
421: }
422: };
424: // An Overlap is a Sifter describing the overlap of two Sieves
425: // Each arrow is local point ---(remote point)---> remote rank right now
426: // For XSifter, this should change to (local patch, local point) ---> (remote rank, remote patch, remote point)
427: }
429: #endif