Actual source code: Distribution.hh
1: #ifndef included_ALE_Distribution_hh
2: #define included_ALE_Distribution_hh
4: #ifndef included_ALE_Mesh_hh
5: #include <Mesh.hh>
6: #endif
8: #ifndef included_ALE_Completion_hh
9: #include <Completion.hh>
10: #endif
12: // Attempt to unify all of the distribution mechanisms:
13: // one to many (distributeMesh)
14: // many to one (unifyMesh)
15: // many to many (Numbering)
16: // as well as things being distributed
17: // Section
18: // Sieve (This sends two sections, the points and cones)
19: // Numbering (Should be an integer section)
20: // Global Order (should be an integer section with extra methods)
21: //
22: // 0) Create the new object to hold the communicated data
23: //
24: // 1) Create Overlap
25: // There may be special ways to do this based upon what we know at the time
26: //
27: // 2) Create send and receive sections over the interface
28: // These have a flat topology now, consisting only of the overlap nodes
29: // We could make a full topology on the overlap (maybe it is necessary for higher order)
30: //
31: // 3) Communication section
32: // Create sizer sections on interface (uses constant sizer)
33: // Communicate sizes on interface (uses custom filler)
34: // Fill send section
35: // sendSection->startCommunication();
36: // recvSection->startCommunication();
37: // sendSection->endCommunication();
38: // recvSection->endCommunication();
39: //
40: // Create section on interface (uses previous sizer)
41: // Communicate values on interface (uses custom filler)
42: // Same stuff as above
43: //
44: // 4) Update new section with old local values (can be done in between the communication?)
45: // Loop over patches in new topology
46: // Loop over chart from patch in old atlas
47: // If this point is in the new sieve from patch
48: // Set to old fiber dimension
49: // Order and allocate new section
50: // Repeat loop, but update values
51: //
52: // 5) Update new section with old received values
53: // Loop over patches in discrete topology of receive section (these are ranks)
54: // Loop over base of discrete sieve (we should transform this to a chart to match above)
55: // Get new patch from overlap, or should the receive patches be <rank, patch>?
56: // Guaranteed to be in the new sieve from patch (but we could check anyway)
57: // Set to recevied fiber dimension
58: // Order and allocate new section
59: // Repeat loop, but update values
60: //
61: // 6) Synchronize PETSc tags (can I get around this?)
62: namespace ALE {
63: template<typename Mesh, typename Partitioner = ALE::Partitioner<> >
64: class DistributionNew {
65: public:
66: typedef Partitioner partitioner_type;
67: typedef typename Mesh::point_type point_type;
68: typedef OrientedPoint<point_type> oriented_point_type;
69: typedef typename Partitioner::part_type rank_type;
70: typedef ALE::ISection<rank_type, point_type> partition_type;
71: typedef ALE::Section<ALE::Pair<int, point_type>, point_type> cones_type;
72: typedef ALE::Section<ALE::Pair<int, point_type>, oriented_point_type> oriented_cones_type;
73: public:
74: template<typename Sieve, typename NewSieve, typename Renumbering, typename SendOverlap, typename RecvOverlap>
75: static Obj<cones_type> completeCones(const Obj<Sieve>& sieve, const Obj<NewSieve>& newSieve, Renumbering& renumbering, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap) {
76: typedef ALE::ConeSection<Sieve> cones_wrapper_type;
77: Obj<cones_wrapper_type> cones = new cones_wrapper_type(sieve);
78: Obj<cones_type> overlapCones = new cones_type(sieve->comm(), sieve->debug());
80: ALE::Pullback::SimpleCopy::copy(sendMeshOverlap, recvMeshOverlap, cones, overlapCones);
81: if (sieve->debug()) {overlapCones->view("Overlap Cones");}
82: // Inserts cones into parallelMesh (must renumber here)
83: ALE::Pullback::InsertionBinaryFusion::fuse(overlapCones, recvMeshOverlap, renumbering, newSieve);
84: return overlapCones;
85: };
86: template<typename Sieve, typename NewSieve, typename SendOverlap, typename RecvOverlap>
87: static Obj<oriented_cones_type> completeConesV(const Obj<Sieve>& sieve, const Obj<NewSieve>& newSieve, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap) {
88: typedef ALE::OrientedConeSectionV<Sieve> oriented_cones_wrapper_type;
89: Obj<oriented_cones_wrapper_type> cones = new oriented_cones_wrapper_type(sieve);
90: Obj<oriented_cones_type> overlapCones = new oriented_cones_type(sieve->comm(), sieve->debug());
92: ALE::Pullback::SimpleCopy::copy(sendMeshOverlap, recvMeshOverlap, cones, overlapCones);
93: if (sieve->debug()) {overlapCones->view("Overlap Cones");}
94: ALE::Pullback::InsertionBinaryFusion::fuse(overlapCones, recvMeshOverlap, newSieve);
95: return overlapCones;
96: };
97: template<typename Sieve, typename NewSieve, typename Renumbering, typename SendOverlap, typename RecvOverlap>
98: static Obj<oriented_cones_type> completeConesV(const Obj<Sieve>& sieve, const Obj<NewSieve>& newSieve, Renumbering& renumbering, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap) {
99: typedef ALE::OrientedConeSectionV<Sieve> oriented_cones_wrapper_type;
100: Obj<oriented_cones_wrapper_type> cones = new oriented_cones_wrapper_type(sieve);
101: Obj<oriented_cones_type> overlapCones = new oriented_cones_type(sieve->comm(), sieve->debug());
103: ALE::Pullback::SimpleCopy::copy(sendMeshOverlap, recvMeshOverlap, cones, overlapCones);
104: if (sieve->debug()) {overlapCones->view("Overlap Cones");}
105: // Inserts cones into parallelMesh (must renumber here)
106: ALE::Pullback::InsertionBinaryFusion::fuse(overlapCones, recvMeshOverlap, renumbering, newSieve);
107: return overlapCones;
108: };
109: // Given a partition of sieve points, copy the mesh pieces to each process and fuse into the new mesh
110: // Return overlaps for the cone communication
111: template<typename Renumbering, typename NewMesh, typename SendOverlap, typename RecvOverlap>
112: static void completeMesh(const Obj<Mesh>& mesh, const Obj<partition_type>& partition, Renumbering& renumbering, const Obj<NewMesh>& newMesh, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap) {
113: typedef ALE::Sifter<rank_type,rank_type,rank_type> part_send_overlap_type;
114: typedef ALE::Sifter<rank_type,rank_type,rank_type> part_recv_overlap_type;
115: const Obj<part_send_overlap_type> sendOverlap = new part_send_overlap_type(partition->comm());
116: const Obj<part_recv_overlap_type> recvOverlap = new part_recv_overlap_type(partition->comm());
118: // Create overlap for partition points
119: // TODO: This needs to be generalized for multiple sources
120: Partitioner::createDistributionPartOverlap(sendOverlap, recvOverlap);
121: // Communicate partition pieces to processes
122: Obj<partition_type> overlapPartition = new partition_type(partition->comm(), partition->debug());
124: overlapPartition->setChart(partition->getChart());
125: ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, partition, overlapPartition);
126: // Create renumbering
127: const int rank = partition->commRank();
128: const point_type *localPoints = partition->restrictPoint(rank);
129: const int numLocalPoints = partition->getFiberDimension(rank);
131: for(point_type p = 0; p < numLocalPoints; ++p) {
132: renumbering[localPoints[p]] = p;
133: }
134: const Obj<typename part_recv_overlap_type::traits::baseSequence> rPoints = recvOverlap->base();
135: point_type localPoint = numLocalPoints;
137: for(typename part_recv_overlap_type::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rPoints->end(); ++p_iter) {
138: const Obj<typename part_recv_overlap_type::coneSequence>& ranks = recvOverlap->cone(*p_iter);
139: const rank_type& remotePartPoint = ranks->begin().color();
140: const point_type *points = overlapPartition->restrictPoint(remotePartPoint);
141: const int numPoints = overlapPartition->getFiberDimension(remotePartPoint);
143: for(int i = 0; i < numPoints; ++i) {
144: renumbering[points[i]] = localPoint++;
145: }
146: }
147: // Create mesh overlap from partition overlap
148: // TODO: Generalize to redistribution (receive from multiple sources)
149: Partitioner::createDistributionMeshOverlap(partition, recvOverlap, renumbering, overlapPartition, sendMeshOverlap, recvMeshOverlap);
150: // Send cones
151: completeCones(mesh->getSieve(), newMesh->getSieve(), renumbering, sendMeshOverlap, recvMeshOverlap);
152: };
153: template<typename Renumbering, typename NewMesh, typename SendOverlap, typename RecvOverlap>
154: static void completeBaseV(const Obj<Mesh>& mesh, const Obj<partition_type>& partition, Renumbering& renumbering, const Obj<NewMesh>& newMesh, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap) {
155: typedef ALE::Sifter<rank_type,rank_type,rank_type> part_send_overlap_type;
156: typedef ALE::Sifter<rank_type,rank_type,rank_type> part_recv_overlap_type;
157: const Obj<part_send_overlap_type> sendOverlap = new part_send_overlap_type(partition->comm());
158: const Obj<part_recv_overlap_type> recvOverlap = new part_recv_overlap_type(partition->comm());
160: // Create overlap for partition points
161: // TODO: This needs to be generalized for multiple sources
162: Partitioner::createDistributionPartOverlap(sendOverlap, recvOverlap);
163: // Communicate partition pieces to processes
164: Obj<partition_type> overlapPartition = new partition_type(partition->comm(), partition->debug());
166: overlapPartition->setChart(partition->getChart());
167: ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, partition, overlapPartition);
168: // Create renumbering
169: const int rank = partition->commRank();
170: const point_type *localPoints = partition->restrictPoint(rank);
171: const int numLocalPoints = partition->getFiberDimension(rank);
173: for(point_type p = 0; p < numLocalPoints; ++p) {
174: ///std::cout <<"["<<partition->commRank()<<"]: local renumbering " << localPoints[p] << " --> " << p << std::endl;
175: renumbering[localPoints[p]] = p;
176: }
177: const Obj<typename part_recv_overlap_type::traits::baseSequence> rPoints = recvOverlap->base();
178: point_type localPoint = numLocalPoints;
180: for(typename part_recv_overlap_type::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rPoints->end(); ++p_iter) {
181: const Obj<typename part_recv_overlap_type::coneSequence>& ranks = recvOverlap->cone(*p_iter);
182: const rank_type& remotePartPoint = ranks->begin().color();
183: const point_type *points = overlapPartition->restrictPoint(remotePartPoint);
184: const int numPoints = overlapPartition->getFiberDimension(remotePartPoint);
186: for(int i = 0; i < numPoints; ++i) {
187: ///std::cout <<"["<<partition->commRank()<<"]: remote renumbering " << points[i] << " --> " << localPoint << std::endl;
188: renumbering[points[i]] = localPoint++;
189: }
190: }
191: newMesh->getSieve()->setChart(typename NewMesh::sieve_type::chart_type(0, renumbering.size()));
192: // Create mesh overlap from partition overlap
193: // TODO: Generalize to redistribution (receive from multiple sources)
194: Partitioner::createDistributionMeshOverlap(partition, recvOverlap, renumbering, overlapPartition, sendMeshOverlap, recvMeshOverlap);
195: };
196: template<typename NewMesh, typename Renumbering, typename SendOverlap, typename RecvOverlap>
197: static Obj<partition_type> distributeMesh(const Obj<Mesh>& mesh, const Obj<NewMesh>& newMesh, Renumbering& renumbering, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap, const int height = 0) {
198: const Obj<partition_type> cellPartition = new partition_type(mesh->comm(), 0, mesh->commSize(), mesh->debug());
199: const Obj<partition_type> partition = new partition_type(mesh->comm(), 0, mesh->commSize(), mesh->debug());
201: // Create the cell partition
202: Partitioner::createPartition(mesh, cellPartition, height);
203: if (mesh->debug()) {
204: PetscViewer viewer;
207: cellPartition->view("Cell Partition");
208: PetscViewerCreate(mesh->comm(), &viewer);CHKERRXX(ierr);
209: PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);CHKERRXX(ierr);
210: PetscViewerFileSetName(viewer, "mesh.vtk");CHKERRXX(ierr);
211: ///TODO MeshView_Sieve_Ascii(mesh, cellPartition, viewer);CHKERRXX(ierr);
212: PetscViewerDestroy(viewer);CHKERRXX(ierr);
213: }
214: // Close the partition over sieve points
215: Partitioner::createPartitionClosure(mesh, cellPartition, partition, height);
216: if (mesh->debug()) {partition->view("Partition");}
217: // Create the remote meshes
218: completeMesh(mesh, partition, renumbering, newMesh, sendMeshOverlap, recvMeshOverlap);
219: // Create the local mesh
220: Partitioner::createLocalMesh(mesh, partition, renumbering, newMesh, height);
221: newMesh->stratify();
222: return partition;
223: };
224: template<typename NewMesh, typename Renumbering, typename SendOverlap, typename RecvOverlap>
225: static Obj<partition_type> distributeMeshAndSections(const Obj<Mesh>& mesh, const Obj<NewMesh>& newMesh, Renumbering& renumbering, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap, const int height = 0) {
226: Obj<partition_type> partition = distributeMesh(mesh, newMesh, renumbering, sendMeshOverlap, recvMeshOverlap, height);
228: // Distribute the coordinates
229: const Obj<typename Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
230: const Obj<typename Mesh::real_section_type>& parallelCoordinates = newMesh->getRealSection("coordinates");
232: newMesh->setupCoordinates(parallelCoordinates);
233: distributeSection(coordinates, partition, renumbering, sendMeshOverlap, recvMeshOverlap, parallelCoordinates);
234: // Distribute other sections
235: if (mesh->getRealSections()->size() > 1) {
236: Obj<std::set<std::string> > names = mesh->getRealSections();
238: for(std::set<std::string>::const_iterator n_iter = names->begin(); n_iter != names->end(); ++n_iter) {
239: if (*n_iter == "coordinates") continue;
240: distributeSection(mesh->getRealSection(*n_iter), partition, renumbering, sendMeshOverlap, recvMeshOverlap, newMesh->getRealSection(*n_iter));
241: }
242: }
243: if (mesh->getIntSections()->size() > 0) {
244: Obj<std::set<std::string> > names = mesh->getIntSections();
246: for(std::set<std::string>::const_iterator n_iter = names->begin(); n_iter != names->end(); ++n_iter) {
247: distributeSection(mesh->getIntSection(*n_iter), partition, renumbering, sendMeshOverlap, recvMeshOverlap, newMesh->getIntSection(*n_iter));
248: }
249: }
250: if (mesh->getArrowSections()->size() > 1) {
251: throw ALE::Exception("Need to distribute more arrow sections");
252: }
253: // Distribute labels
254: const typename Mesh::labels_type& labels = mesh->getLabels();
256: for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
257: if (newMesh->hasLabel(l_iter->first)) continue;
258: const Obj<typename Mesh::label_type>& origLabel = l_iter->second;
259: const Obj<typename Mesh::label_type>& newLabel = newMesh->createLabel(l_iter->first);
260: // Get remote labels
261: ALE::New::Completion<Mesh,typename Mesh::point_type>::scatterCones(origLabel, newLabel, sendMeshOverlap, recvMeshOverlap, renumbering);
262: // Create local label
263: newLabel->add(origLabel, newMesh->getSieve(), renumbering);
264: }
265: return partition;
266: };
267: template<typename NewMesh, typename Renumbering, typename SendOverlap, typename RecvOverlap>
268: static Obj<partition_type> distributeMeshV(const Obj<Mesh>& mesh, const Obj<NewMesh>& newMesh, Renumbering& renumbering, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap, const int height = 0) {
269: const Obj<partition_type> cellPartition = new partition_type(mesh->comm(), 0, mesh->commSize(), mesh->debug());
270: const Obj<partition_type> partition = new partition_type(mesh->comm(), 0, mesh->commSize(), mesh->debug());
272: PETSc::Log::Event("DistributeMesh").begin();
273: // Create the cell partition
274: Partitioner::createPartitionV(mesh, cellPartition, height);
275: if (mesh->debug()) {
276: PetscViewer viewer;
279: cellPartition->view("Cell Partition");
280: PetscViewerCreate(mesh->comm(), &viewer);CHKERRXX(ierr);
281: PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);CHKERRXX(ierr);
282: PetscViewerFileSetName(viewer, "mesh.vtk");CHKERRXX(ierr);
283: ///TODO MeshView_Sieve_Ascii(mesh, cellPartition, viewer);CHKERRXX(ierr);
284: PetscViewerDestroy(viewer);CHKERRXX(ierr);
285: }
286: // Close the partition over sieve points
287: Partitioner::createPartitionClosureV(mesh, cellPartition, partition, height);
288: if (mesh->debug()) {partition->view("Partition");}
289: // Create the remote bases
290: completeBaseV(mesh, partition, renumbering, newMesh, sendMeshOverlap, recvMeshOverlap);
291: // Size the local mesh
292: Partitioner::sizeLocalMeshV(mesh, partition, renumbering, newMesh, height);
293: // Create the remote meshes
294: completeConesV(mesh->getSieve(), newMesh->getSieve(), renumbering, sendMeshOverlap, recvMeshOverlap);
295: // Create the local mesh
296: Partitioner::createLocalMeshV(mesh, partition, renumbering, newMesh, height);
297: newMesh->getSieve()->symmetrize();
298: newMesh->stratify();
299: PETSc::Log::Event("DistributeMesh").end();
300: return partition;
301: };
302: // distributeMeshV:
303: // createPartitionV (can be dumb)
304: // createPartitionClosureV (should be low memory)
305: // completeBaseV (???)
306: // Partitioner::createDistributionPartOverlap (low memory)
307: // copy points to partitions (uses small overlap and fake sections)
308: // renumber (map is potentially big, can measure)
309: // Partitioner::createDistributionMeshOverlap (should be large for distribution)
310: // sendMeshOverlap is localPoint--- remotePoint --->remoteRank
311: // recvMeshOverlap is remoteRank--- remotePoint --->localPoint
312: // sizeLocalMeshV (should be low memory)
313: // completeConesV (???)
314: // createLocalMesh (should be low memory)
315: // symmetrize
316: // stratify
317: template<typename NewMesh>
318: static void distributeMeshAndSectionsV(const Obj<Mesh>& mesh, const Obj<NewMesh>& newMesh) {
319: typedef typename Mesh::point_type point_type;
321: const Obj<typename Mesh::send_overlap_type> sendMeshOverlap = new typename Mesh::send_overlap_type(mesh->comm(), mesh->debug());
322: const Obj<typename Mesh::recv_overlap_type> recvMeshOverlap = new typename Mesh::recv_overlap_type(mesh->comm(), mesh->debug());
323: std::map<point_type,point_type>& renumbering = newMesh->getRenumbering();
324: // Distribute the mesh
325: Obj<partition_type> partition = distributeMeshV(mesh, newMesh, renumbering, sendMeshOverlap, recvMeshOverlap);
326: if (mesh->debug()) {
327: std::cout << "["<<mesh->commRank()<<"]: Mesh Renumbering:" << std::endl;
328: for(typename Mesh::renumbering_type::const_iterator r_iter = renumbering.begin(); r_iter != renumbering.end(); ++r_iter) {
329: std::cout << "["<<mesh->commRank()<<"]: global point " << r_iter->first << " --> " << " local point " << r_iter->second << std::endl;
330: }
331: }
332: // Distribute the coordinates
333: PETSc::Log::Event("DistributeCoords").begin();
334: const Obj<typename Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
335: const Obj<typename Mesh::real_section_type>& parallelCoordinates = newMesh->getRealSection("coordinates");
337: newMesh->setupCoordinates(parallelCoordinates);
338: distributeSection(coordinates, partition, renumbering, sendMeshOverlap, recvMeshOverlap, parallelCoordinates);
339: PETSc::Log::Event("DistributeCoords").end();
340: // Distribute other sections
341: if (mesh->getRealSections()->size() > 1) {
342: PETSc::Log::Event("DistributeRealSec").begin();
343: Obj<std::set<std::string> > names = mesh->getRealSections();
344: int n = 0;
346: for(std::set<std::string>::const_iterator n_iter = names->begin(); n_iter != names->end(); ++n_iter) {
347: if (*n_iter == "coordinates") continue;
348: std::cout << "ERROR: Did not distribute real section " << *n_iter << std::endl;
349: ++n;
350: }
351: PETSc::Log::Event("DistributeRealSec").end();
352: if (n) {throw ALE::Exception("Need to distribute more real sections");}
353: }
354: if (mesh->getIntSections()->size() > 0) {
355: PETSc::Log::Event("DistributeIntSec").begin();
356: Obj<std::set<std::string> > names = mesh->getIntSections();
358: for(std::set<std::string>::const_iterator n_iter = names->begin(); n_iter != names->end(); ++n_iter) {
359: const Obj<typename Mesh::int_section_type>& section = mesh->getIntSection(*n_iter);
360: const Obj<typename Mesh::int_section_type>& newSection = newMesh->getIntSection(*n_iter);
361:
362: // We assume all integer sections are complete sections
363: newSection->setChart(newMesh->getSieve()->getChart());
364: distributeSection(section, partition, renumbering, sendMeshOverlap, recvMeshOverlap, newSection);
365: }
366: PETSc::Log::Event("DistributeIntSec").end();
367: }
368: if (mesh->getArrowSections()->size() > 1) {
369: throw ALE::Exception("Need to distribute more arrow sections");
370: }
371: // Distribute labels
372: PETSc::Log::Event("DistributeLabels").begin();
373: const typename Mesh::labels_type& labels = mesh->getLabels();
375: for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
376: if (newMesh->hasLabel(l_iter->first)) continue;
377: const Obj<typename Mesh::label_type>& origLabel = l_iter->second;
378: const Obj<typename Mesh::label_type>& newLabel = newMesh->createLabel(l_iter->first);
380: #ifdef IMESH_NEW_LABELS
381: newLabel->setChart(newMesh->getSieve()->getChart());
382: // Size the local mesh
383: Partitioner::sizeLocalSieveV(origLabel, partition, renumbering, newLabel);
384: // Create the remote meshes
385: completeConesV(origLabel, newLabel, renumbering, sendMeshOverlap, recvMeshOverlap);
386: // Create the local mesh
387: Partitioner::createLocalSieveV(origLabel, partition, renumbering, newLabel);
388: newLabel->symmetrize();
389: #else
390: distributeLabelV(newMesh->getSieve(), origLabel, partition, renumbering, sendMeshOverlap, recvMeshOverlap, newLabel);
391: #endif
392: }
393: PETSc::Log::Event("DistributeLabels").end();
394: // Create the parallel overlap
395: PETSc::Log::Event("CreateOverlap").begin();
396: Obj<typename Mesh::send_overlap_type> sendParallelMeshOverlap = newMesh->getSendOverlap();
397: Obj<typename Mesh::recv_overlap_type> recvParallelMeshOverlap = newMesh->getRecvOverlap();
398: // Can I figure this out in a nicer way?
399: ALE::SetFromMap<std::map<point_type,point_type> > globalPoints(renumbering);
401: ALE::OverlapBuilder<>::constructOverlap(globalPoints, renumbering, sendParallelMeshOverlap, recvParallelMeshOverlap);
402: newMesh->setCalculatedOverlap(true);
403: PETSc::Log::Event("CreateOverlap").end();
404: };
405: template<typename Label, typename Partition, typename Renumbering, typename SendOverlap, typename RecvOverlap, typename NewLabel>
406: static void distributeLabel(const Obj<typename Mesh::sieve_type>& sieve, const Obj<Label>& l, const Obj<Partition>& partition, Renumbering& renumbering, const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<NewLabel>& newL) {
407: Partitioner::createLocalSifter(l, partition, renumbering, newL);
408: //completeCones(l, newL, renumbering, sendMeshOverlap, recvMeshOverlap);
409: {
410: typedef ALE::UniformSection<point_type, int> cones_type;
411: typedef ALE::LabelSection<typename Mesh::sieve_type, Label> cones_wrapper_type;
412: Obj<cones_wrapper_type> cones = new cones_wrapper_type(sieve, l);
413: Obj<cones_type> overlapCones = new cones_type(l->comm(), l->debug());
415: ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, cones, overlapCones);
416: if (l->debug()) {overlapCones->view("Overlap Label Values");}
417: // Inserts cones into newL (must renumber here)
418: //ALE::Pullback::InsertionBinaryFusion::fuse(overlapCones, recvOverlap, renumbering, newSieve);
419: {
420: typedef typename cones_type::point_type overlap_point_type;
421: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
422: const typename RecvOverlap::traits::baseSequence::iterator rEnd = rPoints->end();
424: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
425: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
426: const typename RecvOverlap::target_type& localPoint = *p_iter;
427: const typename cones_type::point_type& remotePoint = points->begin().color();
428: const overlap_point_type overlapPoint = overlap_point_type(remotePoint.second, remotePoint.first);
429: const int size = overlapCones->getFiberDimension(overlapPoint);
430: const typename cones_type::value_type *values = overlapCones->restrictPoint(overlapPoint);
432: newL->clearCone(localPoint);
433: for(int i = 0; i < size; ++i) {newL->addCone(values[i], localPoint);}
434: }
435: }
436: }
437: };
438: template<typename Label, typename Partition, typename Renumbering, typename SendOverlap, typename RecvOverlap, typename NewLabel>
439: static void distributeLabelV(const Obj<typename Mesh::sieve_type>& sieve, const Obj<Label>& l, const Obj<Partition>& partition, Renumbering& renumbering, const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<NewLabel>& newL) {
440: Partitioner::createLocalSifter(l, partition, renumbering, newL);
441: //completeCones(l, newL, renumbering, sendMeshOverlap, recvMeshOverlap);
442: {
443: typedef typename Label::alloc_type::template rebind<int>::other alloc_type;
444: typedef LabelBaseSectionV<typename Mesh::sieve_type, Label, alloc_type> atlas_type;
445: typedef ALE::UniformSection<ALE::Pair<int, point_type>, int> cones_type;
446: typedef ALE::LabelSection<typename Mesh::sieve_type, Label, alloc_type, atlas_type> cones_wrapper_type;
447: Obj<cones_wrapper_type> cones = new cones_wrapper_type(sieve, l);
448: Obj<cones_type> overlapCones = new cones_type(l->comm(), l->debug());
450: ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, cones, overlapCones);
451: if (l->debug()) {overlapCones->view("Overlap Label Values");}
452: // Inserts cones into newL (must renumber here)
453: //ALE::Pullback::InsertionBinaryFusion::fuse(overlapCones, recvOverlap, renumbering, newSieve);
454: {
455: typedef typename cones_type::point_type overlap_point_type;
456: const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();
458: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rPoints->end(); ++p_iter) {
459: const Obj<typename RecvOverlap::coneSequence>& points = recvOverlap->cone(*p_iter);
460: const typename RecvOverlap::target_type& localPoint = *p_iter;
461: const typename cones_type::point_type& remotePoint = points->begin().color();
462: const overlap_point_type overlapPoint = overlap_point_type(remotePoint.second, remotePoint.first);
463: const int size = overlapCones->getFiberDimension(overlapPoint);
464: const typename cones_type::value_type *values = overlapCones->restrictPoint(overlapPoint);
466: newL->clearCone(localPoint);
467: for(int i = 0; i < size; ++i) {newL->addCone(values[i], localPoint);}
468: }
469: }
470: }
471: };
472: template<typename Section, typename Partition, typename Renumbering, typename SendOverlap, typename RecvOverlap, typename NewSection>
473: static void distributeSection(const Obj<Section>& s, const Obj<Partition>& partition, Renumbering& renumbering, const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<NewSection>& newS) {
474: Partitioner::createLocalSection(s, partition, renumbering, newS);
475: ALE::Completion::completeSection(sendOverlap, recvOverlap, s, newS);
476: };
477: template<typename NewMesh, typename Renumbering, typename SendOverlap, typename RecvOverlap>
478: static Obj<partition_type> unifyMesh(const Obj<Mesh>& mesh, const Obj<NewMesh>& newMesh, Renumbering& renumbering, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap) {
479: const Obj<partition_type> cellPartition = new partition_type(mesh->comm(), 0, mesh->commSize(), mesh->debug());
480: const Obj<partition_type> partition = new partition_type(mesh->comm(), 0, mesh->commSize(), mesh->debug());
481: const Obj<typename Mesh::label_sequence>& cells = mesh->heightStratum(0);
482: const typename Mesh::label_sequence::iterator cEnd = cells->end();
483: typename Mesh::point_type *values = new typename Mesh::point_type[cells->size()];
484: int c = 0;
486: cellPartition->setFiberDimension(0, cells->size());
487: cellPartition->allocatePoint();
488: for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter, ++c) {
489: values[c] = *c_iter;
490: }
491: cellPartition->updatePoint(0, values);
492: delete [] values;
493: // Close the partition over sieve points
494: Partitioner::createPartitionClosure(mesh, cellPartition, partition);
495: // Create the remote meshes
496: completeMesh(mesh, partition, renumbering, newMesh, sendMeshOverlap, recvMeshOverlap);
497: // Create the local mesh
498: Partitioner::createLocalMesh(mesh, partition, renumbering, newMesh);
499: newMesh->stratify();
500: newMesh->view("Unified mesh");
501: return partition;
502: };
503: static Obj<Mesh> unifyMesh(const Obj<Mesh>& mesh) {
504: typedef ALE::Sifter<point_type,rank_type,point_type> mesh_send_overlap_type;
505: typedef ALE::Sifter<rank_type,point_type,point_type> mesh_recv_overlap_type;
506: const Obj<Mesh> newMesh = new Mesh(mesh->comm(), mesh->getDimension(), mesh->debug());
507: const Obj<typename Mesh::sieve_type> newSieve = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
508: const Obj<mesh_send_overlap_type> sendMeshOverlap = new mesh_send_overlap_type(mesh->comm(), mesh->debug());
509: const Obj<mesh_recv_overlap_type> recvMeshOverlap = new mesh_recv_overlap_type(mesh->comm(), mesh->debug());
510: std::map<point_type,point_type> renumbering;
512: newMesh->setSieve(newSieve);
513: const Obj<partition_type> partition = unifyMesh(mesh, newMesh, renumbering, sendMeshOverlap, recvMeshOverlap);
514: // Unify coordinates
515: const Obj<typename Mesh::real_section_type>& coordinates = mesh->getRealSection("coordinates");
516: const Obj<typename Mesh::real_section_type>& newCoordinates = newMesh->getRealSection("coordinates");
518: newMesh->setupCoordinates(newCoordinates);
519: distributeSection(coordinates, partition, renumbering, sendMeshOverlap, recvMeshOverlap, newCoordinates);
520: // Unify labels
521: const typename Mesh::labels_type& labels = mesh->getLabels();
523: for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
524: if (newMesh->hasLabel(l_iter->first)) continue;
525: const Obj<typename Mesh::label_type>& label = l_iter->second;
526: const Obj<typename Mesh::label_type>& newLabel = newMesh->createLabel(l_iter->first);
528: //completeCones(label, newLabel, renumbering, sendMeshOverlap, recvMeshOverlap);
529: {
530: typedef ALE::UniformSection<point_type, int> cones_type;
531: typedef ALE::LabelSection<typename Mesh::sieve_type,typename Mesh::label_type> cones_wrapper_type;
532: Obj<cones_wrapper_type> cones = new cones_wrapper_type(mesh->getSieve(), label);
533: Obj<cones_type> overlapCones = new cones_type(label->comm(), label->debug());
535: ALE::Pullback::SimpleCopy::copy(sendMeshOverlap, recvMeshOverlap, cones, overlapCones);
536: if (label->debug()) {overlapCones->view("Overlap Label Values");}
537: // Inserts cones into parallelMesh (must renumber here)
538: //ALE::Pullback::InsertionBinaryFusion::fuse(overlapCones, recvMeshOverlap, renumbering, newSieve);
539: {
540: const Obj<typename mesh_recv_overlap_type::traits::baseSequence> rPoints = recvMeshOverlap->base();
542: for(typename mesh_recv_overlap_type::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rPoints->end(); ++p_iter) {
543: const Obj<typename mesh_recv_overlap_type::coneSequence>& points = recvMeshOverlap->cone(*p_iter);
544: const typename mesh_recv_overlap_type::target_type& localPoint = *p_iter;
545: const typename cones_type::point_type& remotePoint = points->begin().color();
546: const int size = overlapCones->getFiberDimension(remotePoint);
547: const typename cones_type::value_type *values = overlapCones->restrictPoint(remotePoint);
549: newLabel->clearCone(localPoint);
550: for(int i = 0; i < size; ++i) {newLabel->addCone(values[i], localPoint);}
551: }
552: }
553: }
554: //newLabel->add(label, newSieve);
555: Partitioner::createLocalSifter(label, partition, renumbering, newLabel);
556: newLabel->view(l_iter->first.c_str());
557: }
558: return newMesh;
559: };
560: };
561: template<typename Bundle_>
562: class Distribution {
563: public:
564: typedef Bundle_ bundle_type;
565: typedef typename bundle_type::sieve_type sieve_type;
566: typedef typename bundle_type::point_type point_type;
567: typedef typename bundle_type::alloc_type alloc_type;
568: typedef typename bundle_type::send_overlap_type send_overlap_type;
569: typedef typename bundle_type::recv_overlap_type recv_overlap_type;
570: typedef typename ALE::New::Completion<bundle_type, typename sieve_type::point_type> sieveCompletion;
571: typedef typename ALE::New::SectionCompletion<bundle_type, typename bundle_type::real_section_type::value_type> sectionCompletion;
572: public:
575: static void createPartitionOverlap(const Obj<bundle_type>& bundle, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
576: const Obj<send_overlap_type>& topSendOverlap = bundle->getSendOverlap();
577: const Obj<recv_overlap_type>& topRecvOverlap = bundle->getRecvOverlap();
578: const Obj<typename send_overlap_type::traits::baseSequence> base = topSendOverlap->base();
579: const Obj<typename recv_overlap_type::traits::capSequence> cap = topRecvOverlap->cap();
580: const int rank = bundle->commRank();
582: if (base->empty()) {
583: if (rank == 0) {
584: for(int p = 1; p < bundle->commSize(); p++) {
585: // The arrow is from local partition point p (source) to remote partition point p (color) on rank p (target)
586: sendOverlap->addCone(p, p, p);
587: }
588: }
589: } else {
590: for(typename send_overlap_type::traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
591: const int& p = *b_iter;
592: // The arrow is from local partition point p (source) to remote partition point p (color) on rank p (target)
593: sendOverlap->addCone(p, p, p);
594: }
595: }
596: if (cap->empty()) {
597: if (rank != 0) {
598: // The arrow is from local partition point rank (color) on rank 0 (source) to remote partition point rank (target)
599: recvOverlap->addCone(0, rank, rank);
600: }
601: } else {
602: for(typename recv_overlap_type::traits::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
603: const int& p = *c_iter;
604: // The arrow is from local partition point rank (color) on rank p (source) to remote partition point rank (target)
605: recvOverlap->addCone(p, rank, rank);
606: }
607: }
608: };
611: template<typename Partitioner>
612: static typename Partitioner::part_type *createAssignment(const Obj<bundle_type>& bundle, const int dim, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const int height = 0) {
613: // 1) Form partition point overlap a priori
614: createPartitionOverlap(bundle, sendOverlap, recvOverlap);
615: if (bundle->debug()) {
616: sendOverlap->view("Send overlap for partition");
617: recvOverlap->view("Receive overlap for partition");
618: }
619: // 2) Partition the mesh
620: if (height == 0) {
621: return Partitioner::partitionSieve(bundle, dim);
622: } else if (height == 1) {
623: return Partitioner::partitionSieveByFace(bundle, dim);
624: }
625: throw ALE::Exception("Invalid partition height");
626: };
629: // Partition a bundle on process 0 and scatter to all processes
630: static void scatterBundle(const Obj<bundle_type>& bundle, const int dim, const Obj<bundle_type>& bundleNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const std::string& partitioner, const int height = 0, const Obj<bundle_type>& subBundle = NULL, const Obj<bundle_type>& subBundleNew = NULL) {
631: if (height == 0) {
632: if (partitioner == "chaco") {
633: #ifdef PETSC_HAVE_CHACO
634: typedef typename ALE::New::Chaco::Partitioner<bundle_type> Partitioner;
635: typedef typename ALE::New::Partitioner<bundle_type> GenPartitioner;
636: typedef typename Partitioner::part_type part_type;
638: part_type *assignment = scatterBundle<Partitioner>(bundle, dim, bundleNew, sendOverlap, recvOverlap, height);
639: if (!subBundle.isNull() && !subBundleNew.isNull()) {
640: part_type *subAssignment = GenPartitioner::subordinatePartition(bundle, 1, subBundle, assignment);
641: const Obj<sieve_type>& sieve = subBundle->getSieve();
642: const Obj<sieve_type>& sieveNew = new Mesh::sieve_type(subBundle->comm(), subBundle->debug());
643: const int numCells = subBundle->heightStratum(height)->size();
645: subBundleNew->setSieve(sieveNew);
646: sieveCompletion::scatterSieve(subBundle, sieve, dim, sieveNew, sendOverlap, recvOverlap, height, numCells, subAssignment);
647: subBundleNew->stratify();
648: if (subAssignment != NULL) delete [] subAssignment;
649: }
650: if (assignment != NULL) delete [] assignment;
651: #else
652: throw ALE::Exception("Chaco is not installed. Reconfigure with the flag --download-chaco");
653: #endif
654: } else if (partitioner == "parmetis") {
655: #ifdef PETSC_HAVE_PARMETIS
656: typedef typename ALE::New::ParMetis::Partitioner<bundle_type> Partitioner;
657: typedef typename ALE::New::Partitioner<bundle_type> GenPartitioner;
658: typedef typename Partitioner::part_type part_type;
660: part_type *assignment = scatterBundle<Partitioner>(bundle, dim, bundleNew, sendOverlap, recvOverlap, height);
661: if (!subBundle.isNull() && !subBundleNew.isNull()) {
662: part_type *subAssignment = GenPartitioner::subordinatePartition(bundle, 1, subBundle, assignment);
663: const Obj<sieve_type>& sieve = subBundle->getSieve();
664: const Obj<sieve_type>& sieveNew = new Mesh::sieve_type(subBundle->comm(), subBundle->debug());
665: const int numCells = subBundle->heightStratum(height)->size();
667: subBundleNew->setSieve(sieveNew);
668: sieveCompletion::scatterSieve(subBundle, sieve, dim, sieveNew, sendOverlap, recvOverlap, height, numCells, subAssignment);
669: subBundleNew->stratify();
670: if (subAssignment != NULL) delete [] subAssignment;
671: }
672: if (assignment != NULL) delete [] assignment;
673: #else
674: throw ALE::Exception("ParMetis is not installed. Reconfigure with the flag --download-parmetis");
675: #endif
676: } else {
677: throw ALE::Exception("Unknown partitioner");
678: }
679: } else if (height == 1) {
680: if (partitioner == "zoltan") {
681: #ifdef PETSC_HAVE_ZOLTAN
682: typedef typename ALE::New::Zoltan::Partitioner<bundle_type> Partitioner;
683: typedef typename Partitioner::part_type part_type;
685: part_type *assignment = scatterBundle<Partitioner>(bundle, dim, bundleNew, sendOverlap, recvOverlap, height);
686: if (assignment != NULL) delete [] assignment;
687: #else
688: throw ALE::Exception("Zoltan is not installed. Reconfigure with the flag --download-zoltan");
689: #endif
690: } else if (partitioner == "parmetis") {
691: #ifdef PETSC_HAVE_PARMETIS
692: typedef typename ALE::New::ParMetis::Partitioner<bundle_type> Partitioner;
693: typedef typename Partitioner::part_type part_type;
695: part_type *assignment = scatterBundle<Partitioner>(bundle, dim, bundleNew, sendOverlap, recvOverlap, height);
696: if (assignment != NULL) delete [] assignment;
697: #else
698: throw ALE::Exception("ParMetis is not installed. Reconfigure with the flag --download-parmetis");
699: #endif
700: } else {
701: throw ALE::Exception("Unknown partitioner");
702: }
703: } else {
704: throw ALE::Exception("Invalid partition height");
705: }
706: };
707: template<typename Partitioner>
708: static typename Partitioner::part_type *scatterBundle(const Obj<bundle_type>& bundle, const int dim, const Obj<bundle_type>& bundleNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const int height = 0) {
709: typename Partitioner::part_type *assignment = createAssignment<Partitioner>(bundle, dim, sendOverlap, recvOverlap, height);
710: const Obj<sieve_type>& sieve = bundle->getSieve();
711: const Obj<sieve_type>& sieveNew = bundleNew->getSieve();
712: const int numPoints = bundle->heightStratum(height)->size();
714: sieveCompletion::scatterSieve(bundle, sieve, dim, sieveNew, sendOverlap, recvOverlap, height, numPoints, assignment);
715: bundleNew->stratify();
716: return assignment;
717: };
720: static Obj<ALE::Mesh> distributeMesh(const Obj<ALE::Mesh>& serialMesh, const int height = 0, const std::string& partitioner = "chaco") {
721: MPI_Comm comm = serialMesh->comm();
722: const int dim = serialMesh->getDimension();
723: Obj<ALE::Mesh> parallelMesh = new ALE::Mesh(comm, dim, serialMesh->debug());
724: const Obj<ALE::Mesh::sieve_type>& parallelSieve = new ALE::Mesh::sieve_type(comm, serialMesh->debug());
726: ALE_LOG_EVENT_BEGIN;
727: parallelMesh->setSieve(parallelSieve);
728: if (serialMesh->debug()) {serialMesh->view("Serial mesh");}
730: // Distribute cones
731: Obj<send_overlap_type> sendOverlap = new send_overlap_type(comm, serialMesh->debug());
732: Obj<recv_overlap_type> recvOverlap = new recv_overlap_type(comm, serialMesh->debug());
733: scatterBundle(serialMesh, dim, parallelMesh, sendOverlap, recvOverlap, partitioner, height);
734: parallelMesh->setDistSendOverlap(sendOverlap);
735: parallelMesh->setDistRecvOverlap(recvOverlap);
737: // Distribute labels
738: const typename bundle_type::labels_type& labels = serialMesh->getLabels();
740: for(typename bundle_type::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
741: if (parallelMesh->hasLabel(l_iter->first)) continue;
742: const Obj<typename bundle_type::label_type>& serialLabel = l_iter->second;
743: const Obj<typename bundle_type::label_type>& parallelLabel = parallelMesh->createLabel(l_iter->first);
744: // Create local label
745: #define NEW_LABEL
746: #ifdef NEW_LABEL
747: parallelLabel->add(serialLabel, parallelSieve);
748: #else
749: const Obj<typename bundle_type::label_type::traits::baseSequence>& base = serialLabel->base();
751: for(typename bundle_type::label_type::traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
752: if (parallelSieve->capContains(*b_iter) || parallelSieve->baseContains(*b_iter)) {
753: parallelLabel->addArrow(*serialLabel->cone(*b_iter)->begin(), *b_iter);
754: }
755: }
756: #endif
757: // Get remote labels
758: sieveCompletion::scatterCones(serialLabel, parallelLabel, sendOverlap, recvOverlap);
759: }
761: // Distribute sections
762: Obj<std::set<std::string> > sections = serialMesh->getRealSections();
764: for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
765: parallelMesh->setRealSection(*name, distributeSection(serialMesh->getRealSection(*name), parallelMesh, sendOverlap, recvOverlap));
766: }
767: sections = serialMesh->getIntSections();
768: for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
769: parallelMesh->setIntSection(*name, distributeSection(serialMesh->getIntSection(*name), parallelMesh, sendOverlap, recvOverlap));
770: }
771: sections = serialMesh->getArrowSections();
773: for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
774: parallelMesh->setArrowSection(*name, distributeArrowSection(serialMesh->getArrowSection(*name), serialMesh, parallelMesh, sendOverlap, recvOverlap));
775: }
776: if (parallelMesh->debug()) {parallelMesh->view("Parallel Mesh");}
777: ALE_LOG_EVENT_END;
778: return parallelMesh;
779: };
782: template<typename Section>
783: static void updateSectionLocal(const Obj<Section>& oldSection, const Obj<bundle_type>& newBundle, const Obj<Section>& newSection) {
784: const Obj<typename bundle_type::sieve_type>& newSieve = newBundle->getSieve();
785: const typename Section::atlas_type::chart_type& oldChart = oldSection->getChart();
787: for(typename Section::atlas_type::chart_type::const_iterator c_iter = oldChart.begin(); c_iter != oldChart.end(); ++c_iter) {
788: if (newSieve->capContains(*c_iter) || newSieve->baseContains(*c_iter)) {
789: newSection->setFiberDimension(*c_iter, oldSection->getFiberDimension(*c_iter));
790: }
791: }
792: newBundle->allocate(newSection);
793: const typename Section::atlas_type::chart_type& newChart = newSection->getChart();
795: for(typename Section::atlas_type::chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
796: newSection->updatePointAll(*c_iter, oldSection->restrictPoint(*c_iter));
797: }
798: };
801: template<typename RecvSection, typename Section>
802: static void updateSectionRemote(const Obj<recv_overlap_type>& recvOverlap, const Obj<RecvSection>& recvSection, const Obj<bundle_type>& newBundle, const Obj<Section>& newSection) {
803: Obj<typename recv_overlap_type::traits::baseSequence> recvPoints = recvOverlap->base();
805: for(typename recv_overlap_type::traits::baseSequence::iterator r_iter = recvPoints->begin(); r_iter != recvPoints->end(); ++r_iter) {
806: const Obj<typename recv_overlap_type::traits::coneSequence>& recvPatches = recvOverlap->cone(*r_iter);
807: const typename recv_overlap_type::traits::coneSequence::iterator end = recvPatches->end();
809: for(typename recv_overlap_type::traits::coneSequence::iterator p_iter = recvPatches->begin(); p_iter != end; ++p_iter) {
810: newSection->addPoint(*r_iter, recvSection->getSection(*p_iter)->getFiberDimension(*r_iter));
811: }
812: }
813: newBundle->reallocate(newSection);
814: for(typename recv_overlap_type::traits::baseSequence::iterator r_iter = recvPoints->begin(); r_iter != recvPoints->end(); ++r_iter) {
815: const Obj<typename recv_overlap_type::traits::coneSequence>& recvPatches = recvOverlap->cone(*r_iter);
816: const typename recv_overlap_type::traits::coneSequence::iterator end = recvPatches->end();
818: for(typename recv_overlap_type::traits::coneSequence::iterator p_iter = recvPatches->begin(); p_iter != end; ++p_iter) {
819: if (recvSection->getSection(*p_iter)->getFiberDimension(*r_iter)) {
820: newSection->updatePointAll(*r_iter, recvSection->getSection(*p_iter)->restrictPoint(*r_iter));
821: }
822: }
823: }
824: };
827: template<typename Section>
828: static Obj<Section> distributeSection(const Obj<Section>& serialSection, const Obj<bundle_type>& parallelBundle, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
829: if (serialSection->debug()) {
830: serialSection->view("Serial Section");
831: }
832: typedef typename alloc_type::template rebind<typename Section::value_type>::other value_alloc_type;
833: typedef ALE::Field<send_overlap_type, int, ALE::Section<point_type, typename Section::value_type, value_alloc_type> > send_section_type;
834: typedef ALE::Field<recv_overlap_type, int, ALE::Section<point_type, typename Section::value_type, value_alloc_type> > recv_section_type;
835: typedef ALE::New::SizeSection<Section> SectionSizer;
836: Obj<Section> parallelSection = new Section(serialSection->comm(), serialSection->debug());
837: const Obj<send_section_type> sendSection = new send_section_type(serialSection->comm(), serialSection->debug());
838: const Obj<recv_section_type> recvSection = new recv_section_type(serialSection->comm(), sendSection->getTag(), serialSection->debug());
839: const Obj<SectionSizer> sizer = new SectionSizer(serialSection);
841: updateSectionLocal(serialSection, parallelBundle, parallelSection);
842: sectionCompletion::completeSection(sendOverlap, recvOverlap, sizer, serialSection, sendSection, recvSection);
843: updateSectionRemote(recvOverlap, recvSection, parallelBundle, parallelSection);
844: if (parallelSection->debug()) {
845: parallelSection->view("Parallel Section");
846: }
847: return parallelSection;
848: };
851: template<typename Section>
852: static void updateArrowSectionLocal(const Obj<Section>& oldSection, const Obj<bundle_type>& newBundle, const Obj<Section>& newSection) {
853: const Obj<typename bundle_type::sieve_type>& newSieve = newBundle->getSieve();
854: const typename Section::atlas_type::chart_type& oldChart = oldSection->getChart();
856: for(typename Section::atlas_type::chart_type::const_iterator c_iter = oldChart.begin(); c_iter != oldChart.end(); ++c_iter) {
857: // Dmitry should provide a Sieve::contains(MinimalArrow) method
858: if (newSieve->capContains(c_iter->source) && newSieve->baseContains(c_iter->target)) {
859: newSection->setFiberDimension(*c_iter, oldSection->getFiberDimension(*c_iter));
860: }
861: }
862: //newBundle->allocate(newSection);
863: const typename Section::atlas_type::chart_type& newChart = newSection->getChart();
865: for(typename Section::atlas_type::chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
866: newSection->updatePointAll(*c_iter, oldSection->restrictPoint(*c_iter));
867: }
868: };
871: template<typename RecvSection, typename Section>
872: static void updateArrowSectionRemote(const Obj<recv_overlap_type>& recvOverlap, const Obj<RecvSection>& recvSection, const Obj<bundle_type>& newBundle, const Obj<Section>& newSection) {
873: Obj<typename recv_overlap_type::traits::baseSequence> recvPoints = recvOverlap->base();
875: for(typename recv_overlap_type::traits::baseSequence::iterator r_iter = recvPoints->begin(); r_iter != recvPoints->end(); ++r_iter) {
876: const Obj<typename bundle_type::sieve_type::traits::coneSequence>& cone = newBundle->getSieve()->cone(*r_iter);
877: const typename bundle_type::sieve_type::traits::coneSequence::iterator end = cone->end();
879: for(typename bundle_type::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
880: newSection->setFiberDimension(typename Section::point_type(*c_iter, *r_iter), 1);
881: }
882: }
883: //newBundle->reallocate(newSection);
884: for(typename recv_overlap_type::traits::baseSequence::iterator r_iter = recvPoints->begin(); r_iter != recvPoints->end(); ++r_iter) {
885: const Obj<typename recv_overlap_type::traits::coneSequence>& recvPatches = recvOverlap->cone(*r_iter);
886: const typename recv_overlap_type::traits::coneSequence::iterator recvEnd = recvPatches->end();
888: for(typename recv_overlap_type::traits::coneSequence::iterator p_iter = recvPatches->begin(); p_iter != recvEnd; ++p_iter) {
889: const Obj<typename RecvSection::section_type>& section = recvSection->getSection(*p_iter);
891: if (section->getFiberDimension(*r_iter)) {
892: const Obj<typename bundle_type::sieve_type::traits::coneSequence>& cone = newBundle->getSieve()->cone(*r_iter);
893: const typename bundle_type::sieve_type::traits::coneSequence::iterator end = cone->end();
894: const typename RecvSection::value_type *values = section->restrictPoint(*r_iter);
895: int c = -1;
897: for(typename bundle_type::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
898: newSection->updatePoint(typename Section::point_type(*c_iter, *r_iter), &values[++c]);
899: }
900: }
901: }
902: }
903: };
906: template<typename Section>
907: static Obj<Section> distributeArrowSection(const Obj<Section>& serialSection, const Obj<bundle_type>& serialBundle, const Obj<bundle_type>& parallelBundle, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
908: if (serialSection->debug()) {
909: serialSection->view("Serial ArrowSection");
910: }
911: typedef typename alloc_type::template rebind<typename Section::value_type>::other value_alloc_type;
912: typedef ALE::Field<send_overlap_type, int, ALE::Section<point_type, typename Section::value_type, value_alloc_type> > send_section_type;
913: typedef ALE::Field<recv_overlap_type, int, ALE::Section<point_type, typename Section::value_type, value_alloc_type> > recv_section_type;
914: typedef ALE::New::ConeSizeSection<bundle_type, sieve_type> SectionSizer;
915: typedef ALE::New::ArrowSection<sieve_type, Section> ArrowFiller;
916: Obj<Section> parallelSection = new Section(serialSection->comm(), serialSection->debug());
917: const Obj<send_section_type> sendSection = new send_section_type(serialSection->comm(), serialSection->debug());
918: const Obj<recv_section_type> recvSection = new recv_section_type(serialSection->comm(), sendSection->getTag(), serialSection->debug());
919: const Obj<SectionSizer> sizer = new SectionSizer(serialBundle, serialBundle->getSieve());
920: const Obj<ArrowFiller> filler = new ArrowFiller(serialBundle->getSieve(), serialSection);
922: updateArrowSectionLocal(serialSection, parallelBundle, parallelSection);
923: sectionCompletion::completeSection(sendOverlap, recvOverlap, sizer, filler, sendSection, recvSection);
924: updateArrowSectionRemote(recvOverlap, recvSection, parallelBundle, parallelSection);
925: if (parallelSection->debug()) {
926: parallelSection->view("Parallel ArrowSection");
927: }
928: return parallelSection;
929: };
930: static void unifyBundle(const Obj<bundle_type>& bundle, const int dim, const Obj<bundle_type>& bundleNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
931: typedef int part_type;
932: const Obj<sieve_type>& sieve = bundle->getSieve();
933: const Obj<sieve_type>& sieveNew = bundleNew->getSieve();
934: const int rank = bundle->commRank();
935: const int debug = bundle->debug();
937: // 1) Form partition point overlap a priori
938: if (rank == 0) {
939: for(int p = 1; p < sieve->commSize(); p++) {
940: // The arrow is from remote partition point 0 on rank p to local partition point 0
941: recvOverlap->addCone(p, 0, 0);
942: }
943: } else {
944: // The arrow is from local partition point 0 to remote partition point 0 on rank 0
945: sendOverlap->addCone(0, 0, 0);
946: }
947: if (debug) {
948: sendOverlap->view("Send overlap for partition");
949: recvOverlap->view("Receive overlap for partition");
950: }
951: // 2) Partition the mesh
952: int numCells = bundle->heightStratum(0)->size();
953: part_type *assignment = new part_type[numCells];
955: for(int c = 0; c < numCells; ++c) {
956: assignment[c] = 0;
957: }
958: // 3) Scatter the sieve
959: sieveCompletion::scatterSieve(bundle, sieve, dim, sieveNew, sendOverlap, recvOverlap, 0, numCells, assignment);
960: bundleNew->stratify();
961: // 4) Cleanup
962: if (assignment != NULL) delete [] assignment;
963: };
966: static Obj<ALE::Mesh> unifyMesh(const Obj<ALE::Mesh>& parallelMesh) {
967: const int dim = parallelMesh->getDimension();
968: Obj<ALE::Mesh> serialMesh = new ALE::Mesh(parallelMesh->comm(), dim, parallelMesh->debug());
969: const Obj<ALE::Mesh::sieve_type>& serialSieve = new ALE::Mesh::sieve_type(parallelMesh->comm(), parallelMesh->debug());
971: ALE_LOG_EVENT_BEGIN;
972: serialMesh->setSieve(serialSieve);
973: if (parallelMesh->debug()) {
974: parallelMesh->view("Parallel topology");
975: }
977: // Unify cones
978: Obj<send_overlap_type> sendOverlap = new send_overlap_type(serialMesh->comm(), serialMesh->debug());
979: Obj<recv_overlap_type> recvOverlap = new recv_overlap_type(serialMesh->comm(), serialMesh->debug());
980: unifyBundle(parallelMesh, dim, serialMesh, sendOverlap, recvOverlap);
981: serialMesh->setDistSendOverlap(sendOverlap);
982: serialMesh->setDistRecvOverlap(recvOverlap);
984: // Unify labels
985: const typename bundle_type::labels_type& labels = parallelMesh->getLabels();
987: for(typename bundle_type::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
988: if (serialMesh->hasLabel(l_iter->first)) continue;
989: const Obj<typename bundle_type::label_type>& parallelLabel = l_iter->second;
990: const Obj<typename bundle_type::label_type>& serialLabel = serialMesh->createLabel(l_iter->first);
992: sieveCompletion::scatterCones(parallelLabel, serialLabel, sendOverlap, recvOverlap);
993: }
995: // Unify coordinates
996: Obj<std::set<std::string> > sections = parallelMesh->getRealSections();
998: for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
999: serialMesh->setRealSection(*name, distributeSection(parallelMesh->getRealSection(*name), serialMesh, sendOverlap, recvOverlap));
1000: }
1001: sections = parallelMesh->getIntSections();
1002: for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
1003: serialMesh->setIntSection(*name, distributeSection(parallelMesh->getIntSection(*name), serialMesh, sendOverlap, recvOverlap));
1004: }
1005: sections = parallelMesh->getArrowSections();
1006: for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
1007: serialMesh->setArrowSection(*name, distributeArrowSection(parallelMesh->getArrowSection(*name), parallelMesh, serialMesh, sendOverlap, recvOverlap));
1008: }
1009: if (serialMesh->debug()) {serialMesh->view("Serial Mesh");}
1010: ALE_LOG_EVENT_END;
1011: return serialMesh;
1012: };
1013: public: // Do not like these
1016: // This is just crappy. We could introduce another phase to find out exactly what
1017: // indices people do not have in the global order after communication
1018: template<typename OrigSendOverlap, typename OrigRecvOverlap, typename SendSection, typename RecvSection>
1019: static void updateOverlap(const Obj<OrigSendOverlap>& origSendOverlap, const Obj<OrigRecvOverlap>& origRecvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
1020: const typename SendSection::sheaf_type& sendRanks = sendSection->getPatches();
1021: const typename RecvSection::sheaf_type& recvRanks = recvSection->getPatches();
1023: for(typename SendSection::sheaf_type::const_iterator p_iter = sendRanks.begin(); p_iter != sendRanks.end(); ++p_iter) {
1024: const typename SendSection::patch_type& rank = p_iter->first;
1025: const Obj<typename SendSection::section_type>& section = p_iter->second;
1026: const typename SendSection::section_type::chart_type& chart = section->getChart();
1028: for(typename SendSection::section_type::chart_type::const_iterator b_iter = chart.begin(); b_iter != chart.end(); ++b_iter) {
1029: const typename SendSection::value_type *points = section->restrictPoint(*b_iter);
1030: const int size = section->getFiberDimension(*b_iter);
1032: for(int p = 0; p < size; p++) {
1033: if (origSendOverlap->support(points[p])->size() == 0) {
1034: sendOverlap->addArrow(points[p], rank, points[p]);
1035: }
1036: }
1037: }
1038: }
1039: for(typename RecvSection::sheaf_type::const_iterator p_iter = recvRanks.begin(); p_iter != recvRanks.end(); ++p_iter) {
1040: const typename RecvSection::patch_type& rank = p_iter->first;
1041: const Obj<typename RecvSection::section_type>& section = p_iter->second;
1042: const typename RecvSection::section_type::chart_type& chart = section->getChart();
1044: for(typename RecvSection::section_type::chart_type::const_iterator b_iter = chart.begin(); b_iter != chart.end(); ++b_iter) {
1045: const typename RecvSection::value_type *points = section->restrictPoint(*b_iter);
1046: const int size = section->getFiberDimension(*b_iter);
1048: for(int p = 0; p < size; p++) {
1049: if (origRecvOverlap->support(rank, points[p])->size() == 0) {
1050: recvOverlap->addArrow(rank, points[p], points[p]);
1051: }
1052: }
1053: }
1054: }
1055: };
1058: template<typename RecvOverlap, typename RecvSection>
1059: static void updateSieve(const Obj<RecvOverlap>& recvOverlap, const Obj<RecvSection>& recvSection, const Obj<sieve_type>& sieve) {
1060: #if 1
1061: Obj<typename RecvOverlap::traits::baseSequence> recvPoints = recvOverlap->base();
1063: for(typename RecvOverlap::traits::baseSequence::iterator p_iter = recvPoints->begin(); p_iter != recvPoints->end(); ++p_iter) {
1064: const Obj<typename RecvOverlap::traits::coneSequence>& ranks = recvOverlap->cone(*p_iter);
1065: const typename RecvOverlap::target_type& localPoint = *p_iter;
1067: for(typename RecvOverlap::traits::coneSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
1068: const typename RecvOverlap::target_type& remotePoint = r_iter.color();
1069: const int rank = *r_iter;
1070: const Obj<typename RecvSection::section_type>& section = recvSection->getSection(rank);
1071: const typename RecvSection::value_type *points = section->restrictPoint(remotePoint);
1072: const int size = section->getFiberDimension(remotePoint);
1073: int c = 0;
1075: ///std::cout << "["<<recvSection->commRank()<<"]: Receiving " << size << " points from rank " << rank << std::endl;
1076: for(int p = 0; p < size; p++) {
1077: // rank -- remote point --> local point
1078: if (recvOverlap->support(rank, points[p])->size()) {
1079: sieve->addArrow(*recvOverlap->support(rank, points[p])->begin(), localPoint, c);
1080: ///std::cout << "["<<recvSection->commRank()<<"]: 1Adding arrow " << *recvOverlap->support(rank, points[p])->begin() << "("<<points[p]<<") --> " << localPoint << std::endl;
1081: } else {
1082: sieve->addArrow(points[p], localPoint, c);
1083: ///std::cout << "["<<recvSection->commRank()<<"]: 2Adding arrow " << points[p] << " --> " << localPoint << std::endl;
1084: }
1085: }
1086: }
1087: }
1088: #else
1089: const typename RecvSection::sheaf_type& ranks = recvSection->getPatches();
1091: for(typename RecvSection::sheaf_type::const_iterator p_iter = ranks.begin(); p_iter != ranks.end(); ++p_iter) {
1092: const Obj<typename RecvSection::section_type>& section = p_iter->second;
1093: const typename RecvSection::section_type::chart_type& chart = section->getChart();
1095: for(typename RecvSection::section_type::chart_type::const_iterator b_iter = chart.begin(); b_iter != chart.end(); ++b_iter) {
1096: const typename RecvSection::value_type *points = section->restrictPoint(*b_iter);
1097: int size = section->getFiberDimension(*b_iter);
1098: int c = 0;
1100: std::cout << "["<<recvSection->commRank()<<"]: Receiving " << size << " points from rank " << p_iter->first << std::endl;
1101: for(int p = 0; p < size; p++) {
1102: //sieve->addArrow(points[p], *b_iter, c++);
1103: sieve->addArrow(points[p], *b_iter, c);
1104: std::cout << "["<<recvSection->commRank()<<"]: Adding arrow " << points[p] << " --> " << *b_iter << std::endl;
1105: }
1106: }
1107: }
1108: #endif
1109: };
1112: template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
1113: static void coneCompletion(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<bundle_type>& bundle, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
1114: if (sendOverlap->commSize() == 1) return;
1115: // Distribute cones
1116: const Obj<sieve_type>& sieve = bundle->getSieve();
1117: const Obj<typename sieveCompletion::topology_type> secTopology = sieveCompletion::completion::createSendTopology(sendOverlap);
1118: const Obj<typename sieveCompletion::cone_size_section> coneSizeSection = new typename sieveCompletion::cone_size_section(bundle, sieve);
1119: const Obj<typename sieveCompletion::cone_section> coneSection = new typename sieveCompletion::cone_section(sieve);
1120: sieveCompletion::completion::completeSection(sendOverlap, recvOverlap, coneSizeSection, coneSection, sendSection, recvSection);
1121: // Update cones
1122: updateSieve(recvOverlap, recvSection, sieve);
1123: };
1126: template<typename Section>
1127: static void completeSection(const Obj<bundle_type>& bundle, const Obj<Section>& section) {
1128: typedef typename Distribution<bundle_type>::sieveCompletion sieveCompletion;
1129: typedef typename bundle_type::send_overlap_type send_overlap_type;
1130: typedef typename bundle_type::recv_overlap_type recv_overlap_type;
1131: typedef typename Section::value_type value_type;
1132: typedef typename alloc_type::template rebind<typename Section::value_type>::other value_alloc_type;
1133: typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
1134: typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
1135: typedef ALE::New::SizeSection<Section> SectionSizer;
1136: const int debug = section->debug();
1138: bundle->constructOverlap();
1139: const Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
1140: const Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();
1141: const Obj<send_section_type> sendSection = new send_section_type(section->comm(), section->debug());
1142: const Obj<recv_section_type> recvSection = new recv_section_type(section->comm(), sendSection->getTag(), section->debug());
1143: const Obj<SectionSizer> sizer = new SectionSizer(section);
1145: sectionCompletion::completeSection(sendOverlap, recvOverlap, sizer, section, sendSection, recvSection);
1146: // Update section with remote data
1147: const Obj<typename recv_overlap_type::traits::baseSequence> recvPoints = bundle->getRecvOverlap()->base();
1149: for(typename recv_overlap_type::traits::baseSequence::iterator r_iter = recvPoints->begin(); r_iter != recvPoints->end(); ++r_iter) {
1150: const Obj<typename recv_overlap_type::traits::coneSequence>& recvPatches = recvOverlap->cone(*r_iter);
1151: const typename recv_overlap_type::traits::coneSequence::iterator end = recvPatches->end();
1153: for(typename recv_overlap_type::traits::coneSequence::iterator p_iter = recvPatches->begin(); p_iter != end; ++p_iter) {
1154: if (recvSection->getSection(*p_iter)->getFiberDimension(p_iter.color())) {
1155: if (debug) {std::cout << "["<<section->commRank()<<"]Completed point " << *r_iter << std::endl;}
1156: section->updateAddPoint(*r_iter, recvSection->getSection(*p_iter)->restrictPoint(p_iter.color()));
1157: }
1158: }
1159: }
1160: };
1161: };
1162: }
1163: #endif