Actual source code: Completion.hh
petsc-3.3-p7 2013-05-11
1: #ifndef included_ALE_Completion_hh
2: #define included_ALE_Completion_hh
4: #ifndef included_ALE_Sections_hh
5: #include <sieve/Sections.hh>
6: #endif
8: #ifndef included_ALE_SectionCompletion_hh
9: #include <sieve/SectionCompletion.hh>
10: #endif
12: #ifndef included_ALE_ParallelMapping_hh
13: #include <sieve/ParallelMapping.hh>
14: #endif
16: #include <iostream>
17: #include <fstream>
19: namespace ALE {
20: class Completion {
21: public:
22: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
23: static void completeSection(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
24: typedef ALE::Section<ALE::Pair<int, typename SendOverlap::source_type>, typename SendSection::value_type> OverlapSection;
25: //typedef ALE::Section<typename SendSection::point_type, typename SendSection::value_type> OverlapSection;
26: Obj<OverlapSection> overlapSection = new OverlapSection(sendSection->comm(), sendSection->debug());
28: if (sendSection->debug()) {sendSection->view("Send Section");}
29: ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, sendSection, overlapSection);
30: if (overlapSection->debug()) {overlapSection->view("Overlap Section");}
31: ALE::Pullback::InsertionBinaryFusion::fuse(overlapSection, recvOverlap, recvSection);
32: if (recvSection->debug()) {recvSection->view("Receieve Section");}
33: }
34: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
35: static void completeSectionAdd(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
36: typedef ALE::Section<ALE::Pair<int, typename SendOverlap::source_type>, typename SendSection::value_type> OverlapSection;
37: Obj<OverlapSection> overlapSection = new OverlapSection(sendSection->comm(), sendSection->debug());
39: if (sendSection->debug()) {sendSection->view("Send Section");}
40: ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, sendSection, overlapSection);
41: if (overlapSection->debug()) {overlapSection->view("Overlap Section");}
42: ALE::Pullback::AdditiveBinaryFusion::fuse(overlapSection, recvOverlap, recvSection);
43: if (recvSection->debug()) {recvSection->view("Receieve Section");}
44: }
45: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection, typename OverlapSection>
46: static void completeSectionAdd(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection, const Obj<OverlapSection>& overlapSection) {
47: if (sendSection->debug()) {sendSection->view("Send Section");}
48: ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, sendSection, overlapSection);
49: if (overlapSection->debug()) {overlapSection->view("Overlap Section");}
50: ALE::Pullback::AdditiveBinaryFusion::fuse(overlapSection, recvOverlap, recvSection);
51: if (recvSection->debug()) {recvSection->view("Receieve Section");}
52: }
53: };
54: namespace New {
55: template<typename Bundle_, typename Value_, typename Alloc_ = malloc_allocator<typename Bundle_::point_type> >
56: class Completion {
57: public:
58: typedef int point_type;
59: typedef Value_ value_type;
60: typedef Bundle_ bundle_type;
61: typedef Alloc_ alloc_type;
62: typedef typename alloc_type::template rebind<int>::other int_alloc_type;
63: typedef typename alloc_type::template rebind<value_type>::other value_alloc_type;
64: typedef typename bundle_type::sieve_type sieve_type;
65: typedef typename ALE::DiscreteSieve<point_type, alloc_type> dsieve_type;
66: typedef typename ALE::Topology<int, dsieve_type, alloc_type> topology_type;
67: typedef typename bundle_type::rank_type rank_type;
68: typedef typename ALE::Sifter<point_type, rank_type, point_type> send_overlap_type;
69: typedef typename ALE::Sifter<rank_type, point_type, point_type> recv_overlap_type;
70: typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, int, int_alloc_type> > send_sizer_type;
71: typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, int, int_alloc_type> > recv_sizer_type;
72: typedef typename ALE::New::ConeSizeSection<bundle_type, sieve_type> cone_size_section;
73: typedef typename ALE::New::ConeSection<sieve_type> cone_section;
74: typedef typename ALE::New::SectionCompletion<bundle_type, value_type, alloc_type> completion;
75: public:
76: template<typename PartitionType>
77: static void scatterSieve(const Obj<bundle_type>& bundle, const Obj<sieve_type>& sieve, const int dim, const Obj<sieve_type>& sieveNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const int height, const int numCells, const PartitionType assignment[]) {
78: typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
79: typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
80: int rank = sieve->commRank();
81: int debug = sieve->debug();
83: // Create local sieve
84: const Obj<typename bundle_type::label_sequence>& cells = bundle->heightStratum(height);
85: int e = 0;
87: if (sieve->debug()) {
88: int e2 = 0;
89: for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
90: std::cout << "assignment["<<*e_iter<<"]" << assignment[e2++] << std::endl;
91: }
92: }
93: PetscBool flg;
94: PetscOptionsHasName(PETSC_NULL, "-output_partition", &flg);
95: if (flg) {
96: ostringstream fname;
97: fname << "part." << sieve->commSize() << ".dat";
98: std::ofstream f(fname.str().c_str());
99: int e2 = 0;
100: f << sieve->commSize() << std::endl;
101: for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
102: f << assignment[e2++] << std::endl;
103: }
104: }
105: for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
106: if (assignment[e] == rank) {
107: Obj<typename sieve_type::coneSet> current = new typename sieve_type::coneSet();
108: Obj<typename sieve_type::coneSet> next = new typename sieve_type::coneSet();
109: Obj<typename sieve_type::coneSet> tmp;
111: current->insert(*e_iter);
112: while(current->size()) {
113: for(typename sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
114: const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(*p_iter);
115:
116: for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
117: sieveNew->addArrow(*c_iter, *p_iter, c_iter.color());
118: next->insert(*c_iter);
119: }
120: }
121: tmp = current; current = next; next = tmp;
122: next->clear();
123: }
124: if (height) {
125: current->insert(*e_iter);
126: while(current->size()) {
127: for(typename sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
128: const Obj<typename sieve_type::traits::supportSequence>& support = sieve->support(*p_iter);
129:
130: for(typename sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter) {
131: sieveNew->addArrow(*p_iter, *s_iter, s_iter.color());
132: next->insert(*s_iter);
133: }
134: }
135: tmp = current; current = next; next = tmp;
136: next->clear();
137: }
138: }
139: }
140: e++;
141: }
142: // Complete sizer section
143: typedef typename ALE::New::PartitionSizeSection<bundle_type, PartitionType> partition_size_section;
144: typedef typename ALE::New::PartitionSection<bundle_type, PartitionType> partition_section;
145: Obj<topology_type> secTopology = completion::createSendTopology(sendOverlap);
146: Obj<partition_size_section> partitionSizeSection = new partition_size_section(bundle, height, numCells, assignment);
147: Obj<partition_section> partitionSection = new partition_section(bundle, height, numCells, assignment);
148: Obj<send_section_type> sendSection = new send_section_type(sieve->comm(), sieve->debug());
149: Obj<recv_section_type> recvSection = new recv_section_type(sieve->comm(), sendSection->getTag(), sieve->debug());
151: completion::completeSection(sendOverlap, recvOverlap, partitionSizeSection, partitionSection, sendSection, recvSection);
152: // Unpack the section into the overlap
153: sendOverlap->clear();
154: recvOverlap->clear();
155: const typename send_section_type::sheaf_type& sendPatches = sendSection->getPatches();
157: for(typename send_section_type::sheaf_type::const_iterator p_iter = sendPatches.begin(); p_iter != sendPatches.end(); ++p_iter) {
158: const typename send_section_type::patch_type rank = p_iter->first;
159: const Obj<typename send_section_type::section_type>& section = p_iter->second;
160: const typename send_section_type::section_type::chart_type chart = section->getChart();
162: for(typename send_section_type::section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
163: const typename send_section_type::value_type *points = section->restrictPoint(*c_iter);
164: int size = section->getFiberDimension(*c_iter);
166: for(int p = 0; p < size; p++) {
167: sendOverlap->addArrow(points[p], rank, points[p]);
168: }
169: }
170: }
171: const typename recv_section_type::sheaf_type& recvPatches = recvSection->getPatches();
173: for(typename recv_section_type::sheaf_type::const_iterator p_iter = recvPatches.begin(); p_iter != recvPatches.end(); ++p_iter) {
174: const typename send_section_type::patch_type rank = p_iter->first;
175: const Obj<typename send_section_type::section_type>& section = p_iter->second;
176: const typename send_section_type::section_type::chart_type chart = section->getChart();
178: for(typename recv_section_type::section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
179: const typename recv_section_type::value_type *points = section->restrictPoint(*c_iter);
180: int size = section->getFiberDimension(*c_iter);
182: for(int p = 0; p < size; p++) {
183: recvOverlap->addArrow(rank, points[p], points[p]);
184: }
185: }
186: }
187: if (debug) {
188: sendOverlap->view(std::cout, "Send overlap for points");
189: recvOverlap->view(std::cout, "Receive overlap for points");
190: }
191: // Receive the point section
192: ALE::New::Completion<bundle_type, value_type>::scatterCones(sieve, sieveNew, sendOverlap, recvOverlap, bundle, height);
193: if (height) {
194: ALE::New::Completion<bundle_type, value_type>::scatterSupports(sieve, sieveNew, sendOverlap, recvOverlap, bundle, bundle->depth()-height);
195: }
196: }
197: template<typename SifterType>
198: static void scatterCones(const Obj<SifterType>& sifter, const Obj<SifterType>& sifterNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Obj<bundle_type>& bundle = NULL, const int minimumHeight = 0) {
199: typedef typename ALE::New::ConeSizeSection<bundle_type, SifterType> cone_size_section;
200: typedef typename ALE::New::ConeSection<SifterType> cone_section;
201: typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
202: typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
203: Obj<topology_type> secTopology = completion::createSendTopology(sendOverlap);
204: Obj<cone_size_section> coneSizeSection = new cone_size_section(bundle, sifter, minimumHeight);
205: Obj<cone_section> coneSection = new cone_section(sifter);
206: Obj<send_section_type> sendSection = new send_section_type(sifter->comm(), sifter->debug());
207: Obj<recv_section_type> recvSection = new recv_section_type(sifter->comm(), sendSection->getTag(), sifter->debug());
209: completion::completeSection(sendOverlap, recvOverlap, coneSizeSection, coneSection, sendSection, recvSection);
210: // Unpack the section into the sieve
211: const typename recv_section_type::sheaf_type& patches = recvSection->getPatches();
213: for(typename recv_section_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
214: const Obj<typename recv_section_type::section_type>& section = p_iter->second;
215: const typename recv_section_type::section_type::chart_type& chart = section->getChart();
217: for(typename recv_section_type::section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
218: const typename recv_section_type::value_type *points = section->restrictPoint(*c_iter);
219: int size = section->getFiberDimension(*c_iter);
220: int c = 0;
222: for(int p = 0; p < size; p++) {
223: sifterNew->addArrow(points[p], *c_iter, c++);
224: }
225: }
226: }
227: }
228: template<typename SifterType, typename Renumbering>
229: static void scatterCones(const Obj<SifterType>& sifter, const Obj<SifterType>& sifterNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, Renumbering& renumbering, const Obj<bundle_type>& bundle = NULL) {
230: PETSc::Log::Event("ScatterCones").begin();
231: typedef typename ALE::New::ConeSizeSection<bundle_type, SifterType> cone_size_section;
232: typedef typename ALE::New::ConeSection<SifterType> cone_section;
233: typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
234: typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
235: Obj<topology_type> secTopology = completion::createSendTopology(sendOverlap);
236: Obj<cone_size_section> coneSizeSection = new cone_size_section(bundle, sifter);
237: Obj<cone_section> coneSection = new cone_section(sifter);
238: Obj<send_section_type> sendSection = new send_section_type(sifter->comm(), sifter->debug());
239: Obj<recv_section_type> recvSection = new recv_section_type(sifter->comm(), sendSection->getTag(), sifter->debug());
241: PETSc::Log::Event("ScatterConesComplete").begin();
242: completion::completeSection(sendOverlap, recvOverlap, coneSizeSection, coneSection, sendSection, recvSection);
243: PETSc::Log::Event("ScatterConesComplete").end();
244: PETSc::Log::Event("ScatterConesUpdate").begin();
245: // Unpack the section into the sieve
246: const typename recv_section_type::sheaf_type& patches = recvSection->getPatches();
248: for(typename recv_section_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
249: const Obj<typename recv_section_type::section_type>& section = p_iter->second;
250: const typename recv_section_type::section_type::chart_type& chart = section->getChart();
252: for(typename recv_section_type::section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
253: const typename recv_section_type::value_type *points = section->restrictPoint(*c_iter);
254: int size = section->getFiberDimension(*c_iter);
255: int c = 0;
257: for(int p = 0; p < size; p++) {
258: sifterNew->addArrow(points[p], renumbering[*c_iter], c++);
259: }
260: }
261: }
262: PETSc::Log::Event("ScatterConesUpdate").end();
263: PETSc::Log::Event("ScatterCones").end();
264: }
265: template<typename SifterType>
266: static void scatterSupports(const Obj<SifterType>& sifter, const Obj<SifterType>& sifterNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Obj<bundle_type>& bundle = NULL, const int minimumDepth = 0) {
267: typedef typename ALE::New::SupportSizeSection<bundle_type, SifterType> support_size_section;
268: typedef typename ALE::New::SupportSection<SifterType> support_section;
269: typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
270: typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
271: Obj<topology_type> secTopology = completion::createSendTopology(sendOverlap);
272: Obj<support_size_section> supportSizeSection = new support_size_section(bundle, sifter, minimumDepth);
273: Obj<support_section> supportSection = new support_section(sifter);
274: Obj<send_section_type> sendSection = new send_section_type(sifter->comm(), sifter->debug());
275: Obj<recv_section_type> recvSection = new recv_section_type(sifter->comm(), sendSection->getTag(), sifter->debug());
277: completion::completeSection(sendOverlap, recvOverlap, supportSizeSection, supportSection, sendSection, recvSection);
278: // Unpack the section into the sieve
279: const typename recv_section_type::sheaf_type& recvPatches = recvSection->getPatches();
281: for(typename recv_section_type::sheaf_type::const_iterator p_iter = recvPatches.begin(); p_iter != recvPatches.end(); ++p_iter) {
282: const Obj<typename send_section_type::section_type>& section = p_iter->second;
283: const typename send_section_type::section_type::chart_type chart = section->getChart();
285: for(typename recv_section_type::section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
286: const typename recv_section_type::value_type *points = section->restrictPoint(*c_iter);
287: int size = section->getFiberDimension(*c_iter);
288: int c = 0;
290: for(int p = 0; p < size; p++) {
291: sifterNew->addArrow(*c_iter, points[p], c++);
292: }
293: }
294: }
295: }
296: };
298: template<typename Value_>
299: class ParallelFactory {
300: public:
301: typedef Value_ value_type;
302: protected:
303: int _debug;
304: MPI_Datatype _mpiType;
305: protected:
306: MPI_Datatype constructMPIType() {
307: if (sizeof(value_type) == 1) {
308: return MPI_BYTE;
309: } else if (sizeof(value_type) == 2) {
310: return MPI_SHORT;
311: } else if (sizeof(value_type) == 4) {
312: return MPI_INT;
313: } else if (sizeof(value_type) == 8) {
314: return MPI_DOUBLE;
315: } else if (sizeof(value_type) == 28) {
316: int blen[2];
317: MPI_Aint indices[2];
318: MPI_Datatype oldtypes[2], newtype;
319: blen[0] = 1; indices[0] = 0; oldtypes[0] = MPI_INT;
320: blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
321: MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
322: MPI_Type_commit(&newtype);
323: return newtype;
324: } else if (sizeof(value_type) == 32) {
325: int blen[2];
326: MPI_Aint indices[2];
327: MPI_Datatype oldtypes[2], newtype;
328: blen[0] = 1; indices[0] = 0; oldtypes[0] = MPI_DOUBLE;
329: blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
330: MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
331: MPI_Type_commit(&newtype);
332: return newtype;
333: }
334: throw ALE::Exception("Cannot determine MPI type for value type");
335: };
336: ParallelFactory(const int debug) : _debug(debug) {
337: this->_mpiType = this->constructMPIType();
338: };
339: public:
340: ~ParallelFactory() {};
341: public:
342: static const Obj<ParallelFactory>& singleton(const int debug, bool cleanup = false) {
343: static Obj<ParallelFactory> *_singleton = NULL;
345: if (cleanup) {
346: if (debug) {std::cout << "Destroying ParallelFactory" << std::endl;}
347: if (_singleton) {delete _singleton;}
348: _singleton = NULL;
349: } else if (_singleton == NULL) {
350: if (debug) {std::cout << "Creating new ParallelFactory" << std::endl;}
351: _singleton = new Obj<ParallelFactory>();
352: *_singleton = new ParallelFactory(debug);
353: }
354: return *_singleton;
355: };
356: public: // Accessors
357: int debug() const {return this->_debug;};
358: MPI_Datatype getMPIType() const {return this->_mpiType;};
359: };
360: }
361: }
362: #endif