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