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:       // Create the cell partition
273:       Partitioner::createPartitionV(mesh, cellPartition, height);
274:       if (mesh->debug()) {
275:         PetscViewer    viewer;

278:         cellPartition->view("Cell Partition");
279:         PetscViewerCreate(mesh->comm(), &viewer);CHKERRXX(ierr);
280:         PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);CHKERRXX(ierr);
281:         PetscViewerFileSetName(viewer, "mesh.vtk");CHKERRXX(ierr);
282:         ///TODO MeshView_Sieve_Ascii(mesh, cellPartition, viewer);CHKERRXX(ierr);
283:         PetscViewerDestroy(viewer);CHKERRXX(ierr);
284:       }
285:       // Close the partition over sieve points
286:       Partitioner::createPartitionClosureV(mesh, cellPartition, partition, height);
287:       if (mesh->debug()) {partition->view("Partition");}
288:       // Create the remote bases
289:       completeBaseV(mesh, partition, renumbering, newMesh, sendMeshOverlap, recvMeshOverlap);
290:       // Size the local mesh
291:       Partitioner::sizeLocalMeshV(mesh, partition, renumbering, newMesh, height);
292:       // Create the remote meshes
293:       completeConesV(mesh->getSieve(), newMesh->getSieve(), renumbering, sendMeshOverlap, recvMeshOverlap);
294:       // Create the local mesh
295:       Partitioner::createLocalMeshV(mesh, partition, renumbering, newMesh, height);
296:       newMesh->getSieve()->symmetrize();
297:       newMesh->stratify();
298:       return partition;
299:     };
300:     template<typename NewMesh>
301:     static void distributeMeshAndSectionsV(const Obj<Mesh>& mesh, const Obj<NewMesh>& newMesh) {
302:       typedef typename Mesh::point_type point_type;

304:       const Obj<typename Mesh::send_overlap_type> sendMeshOverlap = new typename Mesh::send_overlap_type(mesh->comm(), mesh->debug());
305:       const Obj<typename Mesh::recv_overlap_type> recvMeshOverlap = new typename Mesh::recv_overlap_type(mesh->comm(), mesh->debug());
306:       std::map<point_type,point_type>&            renumbering     = newMesh->getRenumbering();
307:       // Distribute the mesh
308:       Obj<partition_type> partition = distributeMeshV(mesh, newMesh, renumbering, sendMeshOverlap, recvMeshOverlap);
309:       if (mesh->debug()) {
310:         std::cout << "["<<mesh->commRank()<<"]: Mesh Renumbering:" << std::endl;
311:         for(typename Mesh::renumbering_type::const_iterator r_iter = renumbering.begin(); r_iter != renumbering.end(); ++r_iter) {
312:           std::cout << "["<<mesh->commRank()<<"]:   global point " << r_iter->first << " --> " << " local point " << r_iter->second << std::endl;
313:         }
314:       }
315:       // Distribute the coordinates
316:       const Obj<typename Mesh::real_section_type>& coordinates         = mesh->getRealSection("coordinates");
317:       const Obj<typename Mesh::real_section_type>& parallelCoordinates = newMesh->getRealSection("coordinates");

319:       newMesh->setupCoordinates(parallelCoordinates);
320:       distributeSection(coordinates, partition, renumbering, sendMeshOverlap, recvMeshOverlap, parallelCoordinates);
321:       // Distribute other sections
322:       if (mesh->getRealSections()->size() > 1) {
323:         Obj<std::set<std::string> > names = mesh->getRealSections();
324:         int                         n     = 0;

326:         for(std::set<std::string>::const_iterator n_iter = names->begin(); n_iter != names->end(); ++n_iter) {
327:           if (*n_iter == "coordinates")   continue;
328:           std::cout << "ERROR: Did not distribute real section " << *n_iter << std::endl;
329:           ++n;
330:         }
331:         if (n) {throw ALE::Exception("Need to distribute more real sections");}
332:       }
333:       if (mesh->getIntSections()->size() > 0) {
334:         Obj<std::set<std::string> > names = mesh->getIntSections();

336:         for(std::set<std::string>::const_iterator n_iter = names->begin(); n_iter != names->end(); ++n_iter) {
337:           const Obj<typename Mesh::int_section_type>& section    = mesh->getIntSection(*n_iter);
338:           const Obj<typename Mesh::int_section_type>& newSection = newMesh->getIntSection(*n_iter);
339: 
340:           // We assume all integer sections are complete sections
341:           newSection->setChart(newMesh->getSieve()->getChart());
342:           distributeSection(section, partition, renumbering, sendMeshOverlap, recvMeshOverlap, newSection);
343:         }
344:       }
345:       if (mesh->getArrowSections()->size() > 1) {
346:         throw ALE::Exception("Need to distribute more arrow sections");
347:       }
348:       // Distribute labels
349:       const typename Mesh::labels_type& labels = mesh->getLabels();

351:       for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
352:         if (newMesh->hasLabel(l_iter->first)) continue;
353:         const Obj<typename Mesh::label_type>& origLabel = l_iter->second;
354:         const Obj<typename Mesh::label_type>& newLabel  = newMesh->createLabel(l_iter->first);

356: #ifdef IMESH_NEW_LABELS
357:         newLabel->setChart(newMesh->getSieve()->getChart());
358:         // Size the local mesh
359:         Partitioner::sizeLocalSieveV(origLabel, partition, renumbering, newLabel);
360:         // Create the remote meshes
361:         completeConesV(origLabel, newLabel, renumbering, sendMeshOverlap, recvMeshOverlap);
362:         // Create the local mesh
363:         Partitioner::createLocalSieveV(origLabel, partition, renumbering, newLabel);
364:         newLabel->symmetrize();
365: #else
366:         // Get remote labels
367:         ALE::New::Completion<Mesh,point_type>::scatterCones(origLabel, newLabel, sendMeshOverlap, recvMeshOverlap, renumbering);
368:         // Create local label
369:         newLabel->add(origLabel, newMesh->getSieve(), renumbering);
370: #endif
371:       }
372:       // Create the parallel overlap
373:       Obj<typename Mesh::send_overlap_type> sendParallelMeshOverlap = newMesh->getSendOverlap();
374:       Obj<typename Mesh::recv_overlap_type> recvParallelMeshOverlap = newMesh->getRecvOverlap();
375:       //   Can I figure this out in a nicer way?
376:       ALE::SetFromMap<std::map<point_type,point_type> > globalPoints(renumbering);

378:       ALE::OverlapBuilder<>::constructOverlap(globalPoints, renumbering, sendParallelMeshOverlap, recvParallelMeshOverlap);
379:       newMesh->setCalculatedOverlap(true);
380:     };
381:     template<typename Label, typename Partition, typename Renumbering, typename SendOverlap, typename RecvOverlap, typename NewLabel>
382:     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) {
383:       Partitioner::createLocalSifter(l, partition, renumbering, newL);
384:       //completeCones(l, newL, renumbering, sendMeshOverlap, recvMeshOverlap);
385:       {
386:         typedef ALE::UniformSection<point_type, int>                cones_type;
387:         typedef ALE::LabelSection<typename Mesh::sieve_type, Label> cones_wrapper_type;
388:         Obj<cones_wrapper_type> cones        = new cones_wrapper_type(sieve, l);
389:         Obj<cones_type>         overlapCones = new cones_type(l->comm(), l->debug());

391:         ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, cones, overlapCones);
392:         if (l->debug()) {overlapCones->view("Overlap Label Values");}
393:         // Inserts cones into newL (must renumber here)
394:         //ALE::Pullback::InsertionBinaryFusion::fuse(overlapCones, recvOverlap, renumbering, newSieve);
395:         {
396:           const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();

398:           for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rPoints->end(); ++p_iter) {
399:             const Obj<typename RecvOverlap::coneSequence>& points      = recvOverlap->cone(*p_iter);
400:             const typename RecvOverlap::target_type&       localPoint  = *p_iter;
401:             const typename cones_type::point_type&         remotePoint = points->begin().color();
402:             const int                                      size        = overlapCones->getFiberDimension(remotePoint);
403:             const typename cones_type::value_type         *values      = overlapCones->restrictPoint(remotePoint);

405:             newL->clearCone(localPoint);
406:             for(int i = 0; i < size; ++i) {newL->addCone(values[i], localPoint);}
407:           }
408:         }
409:       }
410:     };
411:     template<typename Section, typename Partition, typename Renumbering, typename SendOverlap, typename RecvOverlap, typename NewSection>
412:     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) {
413:       Partitioner::createLocalSection(s, partition, renumbering, newS);
414:       ALE::Completion::completeSection(sendOverlap, recvOverlap, s, newS);
415:     };
416:     template<typename NewMesh, typename Renumbering, typename SendOverlap, typename RecvOverlap>
417:     static Obj<partition_type> unifyMesh(const Obj<Mesh>& mesh, const Obj<NewMesh>& newMesh, Renumbering& renumbering, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap) {
418:       const Obj<partition_type> cellPartition = new partition_type(mesh->comm(), 0, mesh->commSize(), mesh->debug());
419:       const Obj<partition_type> partition     = new partition_type(mesh->comm(), 0, mesh->commSize(), mesh->debug());
420:       const Obj<typename Mesh::label_sequence>&     cells  = mesh->heightStratum(0);
421:       const typename Mesh::label_sequence::iterator cEnd   = cells->end();
422:       typename Mesh::point_type                    *values = new typename Mesh::point_type[cells->size()];
423:       int                                           c      = 0;

425:       cellPartition->setFiberDimension(0, cells->size());
426:       cellPartition->allocatePoint();
427:       for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter, ++c) {
428:         values[c] = *c_iter;
429:       }
430:       cellPartition->updatePoint(0, values);
431:       delete [] values;
432:       // Close the partition over sieve points
433:       Partitioner::createPartitionClosure(mesh, cellPartition, partition);
434:       // Create the remote meshes
435:       completeMesh(mesh, partition, renumbering, newMesh, sendMeshOverlap, recvMeshOverlap);
436:       // Create the local mesh
437:       Partitioner::createLocalMesh(mesh, partition, renumbering, newMesh);
438:       newMesh->stratify();
439:       newMesh->view("Unified mesh");
440:       return partition;
441:     };
442:     static Obj<Mesh> unifyMesh(const Obj<Mesh>& mesh) {
443:       typedef ALE::Sifter<point_type,rank_type,point_type> mesh_send_overlap_type;
444:       typedef ALE::Sifter<rank_type,point_type,point_type> mesh_recv_overlap_type;
445:       const Obj<Mesh>                      newMesh         = new Mesh(mesh->comm(), mesh->getDimension(), mesh->debug());
446:       const Obj<typename Mesh::sieve_type> newSieve        = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
447:       const Obj<mesh_send_overlap_type>    sendMeshOverlap = new mesh_send_overlap_type(mesh->comm(), mesh->debug());
448:       const Obj<mesh_recv_overlap_type>    recvMeshOverlap = new mesh_recv_overlap_type(mesh->comm(), mesh->debug());
449:       std::map<point_type,point_type>      renumbering;

451:       newMesh->setSieve(newSieve);
452:       const Obj<partition_type> partition = unifyMesh(mesh, newMesh, renumbering, sendMeshOverlap, recvMeshOverlap);
453:       // Unify coordinates
454:       const Obj<typename Mesh::real_section_type>& coordinates    = mesh->getRealSection("coordinates");
455:       const Obj<typename Mesh::real_section_type>& newCoordinates = newMesh->getRealSection("coordinates");

457:       newMesh->setupCoordinates(newCoordinates);
458:       distributeSection(coordinates, partition, renumbering, sendMeshOverlap, recvMeshOverlap, newCoordinates);
459:       // Unify labels
460:       const typename Mesh::labels_type& labels = mesh->getLabels();

462:       for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
463:         if (newMesh->hasLabel(l_iter->first)) continue;
464:         const Obj<typename Mesh::label_type>& label    = l_iter->second;
465:         const Obj<typename Mesh::label_type>& newLabel = newMesh->createLabel(l_iter->first);

467:         //completeCones(label, newLabel, renumbering, sendMeshOverlap, recvMeshOverlap);
468:         {
469:           typedef ALE::UniformSection<point_type, int> cones_type;
470:           typedef ALE::LabelSection<typename Mesh::sieve_type,typename Mesh::label_type> cones_wrapper_type;
471:           Obj<cones_wrapper_type> cones        = new cones_wrapper_type(mesh->getSieve(), label);
472:           Obj<cones_type>         overlapCones = new cones_type(label->comm(), label->debug());

474:           ALE::Pullback::SimpleCopy::copy(sendMeshOverlap, recvMeshOverlap, cones, overlapCones);
475:           if (label->debug()) {overlapCones->view("Overlap Label Values");}
476:           // Inserts cones into parallelMesh (must renumber here)
477:           //ALE::Pullback::InsertionBinaryFusion::fuse(overlapCones, recvMeshOverlap, renumbering, newSieve);
478:           {
479:             const Obj<typename mesh_recv_overlap_type::traits::baseSequence> rPoints = recvMeshOverlap->base();

481:             for(typename mesh_recv_overlap_type::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rPoints->end(); ++p_iter) {
482:               const Obj<typename mesh_recv_overlap_type::coneSequence>& points      = recvMeshOverlap->cone(*p_iter);
483:               const typename mesh_recv_overlap_type::target_type&       localPoint  = *p_iter;
484:               const typename cones_type::point_type&                    remotePoint = points->begin().color();
485:               const int                                                 size        = overlapCones->getFiberDimension(remotePoint);
486:               const typename cones_type::value_type                    *values      = overlapCones->restrictPoint(remotePoint);

488:               newLabel->clearCone(localPoint);
489:               for(int i = 0; i < size; ++i) {newLabel->addCone(values[i], localPoint);}
490:             }
491:           }
492:         }
493:         //newLabel->add(label, newSieve);
494:         Partitioner::createLocalSifter(label, partition, renumbering, newLabel);
495:         newLabel->view(l_iter->first.c_str());
496:       }
497:       return newMesh;
498:     };
499:   };
500:   template<typename Bundle_>
501:   class Distribution {
502:   public:
503:     typedef Bundle_                                                                     bundle_type;
504:     typedef typename bundle_type::sieve_type                                            sieve_type;
505:     typedef typename bundle_type::point_type                                            point_type;
506:     typedef typename bundle_type::alloc_type                                            alloc_type;
507:     typedef typename bundle_type::send_overlap_type                                     send_overlap_type;
508:     typedef typename bundle_type::recv_overlap_type                                     recv_overlap_type;
509:     typedef typename ALE::New::Completion<bundle_type, typename sieve_type::point_type>                            sieveCompletion;
510:     typedef typename ALE::New::SectionCompletion<bundle_type, typename bundle_type::real_section_type::value_type> sectionCompletion;
511:   public:
514:     static void createPartitionOverlap(const Obj<bundle_type>& bundle, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
515:       const Obj<send_overlap_type>& topSendOverlap = bundle->getSendOverlap();
516:       const Obj<recv_overlap_type>& topRecvOverlap = bundle->getRecvOverlap();
517:       const Obj<typename send_overlap_type::traits::baseSequence> base = topSendOverlap->base();
518:       const Obj<typename recv_overlap_type::traits::capSequence>  cap  = topRecvOverlap->cap();
519:       const int rank = bundle->commRank();

521:       if (base->empty()) {
522:         if (rank == 0) {
523:           for(int p = 1; p < bundle->commSize(); p++) {
524:             // The arrow is from local partition point p (source) to remote partition point p (color) on rank p (target)
525:             sendOverlap->addCone(p, p, p);
526:           }
527:         }
528:       } else {
529:         for(typename send_overlap_type::traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
530:           const int& p = *b_iter;
531:           // The arrow is from local partition point p (source) to remote partition point p (color) on rank p (target)
532:           sendOverlap->addCone(p, p, p);
533:         }
534:       }
535:       if (cap->empty()) {
536:         if (rank != 0) {
537:           // The arrow is from local partition point rank (color) on rank 0 (source) to remote partition point rank (target)
538:           recvOverlap->addCone(0, rank, rank);
539:         }
540:       } else {
541:         for(typename recv_overlap_type::traits::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
542:           const int& p = *c_iter;
543:           // The arrow is from local partition point rank (color) on rank p (source) to remote partition point rank (target)
544:           recvOverlap->addCone(p, rank, rank);
545:         }
546:       }
547:     };
550:     template<typename Partitioner>
551:     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) {
552:       // 1) Form partition point overlap a priori
553:       createPartitionOverlap(bundle, sendOverlap, recvOverlap);
554:       if (bundle->debug()) {
555:         sendOverlap->view("Send overlap for partition");
556:         recvOverlap->view("Receive overlap for partition");
557:       }
558:       // 2) Partition the mesh
559:       if (height == 0) {
560:         return Partitioner::partitionSieve(bundle, dim);
561:       } else if (height == 1) {
562:         return Partitioner::partitionSieveByFace(bundle, dim);
563:       }
564:       throw ALE::Exception("Invalid partition height");
565:     };
568:     // Partition a bundle on process 0 and scatter to all processes
569:     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) {
570:       if (height == 0) {
571:         if (partitioner == "chaco") {
572: #ifdef PETSC_HAVE_CHACO
573:           typedef typename ALE::New::Chaco::Partitioner<bundle_type> Partitioner;
574:           typedef typename ALE::New::Partitioner<bundle_type>        GenPartitioner;
575:           typedef typename Partitioner::part_type                    part_type;

577:           part_type *assignment = scatterBundle<Partitioner>(bundle, dim, bundleNew, sendOverlap, recvOverlap, height);
578:           if (!subBundle.isNull() && !subBundleNew.isNull()) {
579:             part_type *subAssignment = GenPartitioner::subordinatePartition(bundle, 1, subBundle, assignment);
580:             const Obj<sieve_type>& sieve      = subBundle->getSieve();
581:             const Obj<sieve_type>& sieveNew   = new Mesh::sieve_type(subBundle->comm(), subBundle->debug());
582:             const int              numCells   = subBundle->heightStratum(height)->size();

584:             subBundleNew->setSieve(sieveNew);
585:             sieveCompletion::scatterSieve(subBundle, sieve, dim, sieveNew, sendOverlap, recvOverlap, height, numCells, subAssignment);
586:             subBundleNew->stratify();
587:             if (subAssignment != NULL) delete [] subAssignment;
588:           }
589:           if (assignment != NULL) delete [] assignment;
590: #else
591:           throw ALE::Exception("Chaco is not installed. Reconfigure with the flag --download-chaco");
592: #endif
593:         } else if (partitioner == "parmetis") {
594: #ifdef PETSC_HAVE_PARMETIS
595:           typedef typename ALE::New::ParMetis::Partitioner<bundle_type> Partitioner;
596:           typedef typename ALE::New::Partitioner<bundle_type>           GenPartitioner;
597:           typedef typename Partitioner::part_type                       part_type;

599:           part_type *assignment = scatterBundle<Partitioner>(bundle, dim, bundleNew, sendOverlap, recvOverlap, height);
600:           if (!subBundle.isNull() && !subBundleNew.isNull()) {
601:             part_type *subAssignment = GenPartitioner::subordinatePartition(bundle, 1, subBundle, assignment);
602:             const Obj<sieve_type>& sieve      = subBundle->getSieve();
603:             const Obj<sieve_type>& sieveNew   = new Mesh::sieve_type(subBundle->comm(), subBundle->debug());
604:             const int              numCells   = subBundle->heightStratum(height)->size();

606:             subBundleNew->setSieve(sieveNew);
607:             sieveCompletion::scatterSieve(subBundle, sieve, dim, sieveNew, sendOverlap, recvOverlap, height, numCells, subAssignment);
608:             subBundleNew->stratify();
609:             if (subAssignment != NULL) delete [] subAssignment;
610:           }
611:           if (assignment != NULL) delete [] assignment;
612: #else
613:           throw ALE::Exception("ParMetis is not installed. Reconfigure with the flag --download-parmetis");
614: #endif
615:         } else {
616:           throw ALE::Exception("Unknown partitioner");
617:         }
618:       } else if (height == 1) {
619:         if (partitioner == "zoltan") {
620: #ifdef PETSC_HAVE_ZOLTAN
621:           typedef typename ALE::New::Zoltan::Partitioner<bundle_type> Partitioner;
622:           typedef typename Partitioner::part_type                     part_type;

624:           part_type *assignment = scatterBundle<Partitioner>(bundle, dim, bundleNew, sendOverlap, recvOverlap, height);
625:           if (assignment != NULL) delete [] assignment;
626: #else
627:           throw ALE::Exception("Zoltan is not installed. Reconfigure with the flag --download-zoltan");
628: #endif
629:         } else if (partitioner == "parmetis") {
630: #ifdef PETSC_HAVE_PARMETIS
631:           typedef typename ALE::New::ParMetis::Partitioner<bundle_type> Partitioner;
632:           typedef typename Partitioner::part_type                       part_type;

634:           part_type *assignment = scatterBundle<Partitioner>(bundle, dim, bundleNew, sendOverlap, recvOverlap, height);
635:           if (assignment != NULL) delete [] assignment;
636: #else
637:           throw ALE::Exception("ParMetis is not installed. Reconfigure with the flag --download-parmetis");
638: #endif
639:         } else {
640:           throw ALE::Exception("Unknown partitioner");
641:         }
642:       } else {
643:         throw ALE::Exception("Invalid partition height");
644:       }
645:     };
646:     template<typename Partitioner>
647:     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) {
648:       typename Partitioner::part_type *assignment = createAssignment<Partitioner>(bundle, dim, sendOverlap, recvOverlap, height);
649:       const Obj<sieve_type>&           sieve      = bundle->getSieve();
650:       const Obj<sieve_type>&           sieveNew   = bundleNew->getSieve();
651:       const int                        numPoints  = bundle->heightStratum(height)->size();

653:       sieveCompletion::scatterSieve(bundle, sieve, dim, sieveNew, sendOverlap, recvOverlap, height, numPoints, assignment);
654:       bundleNew->stratify();
655:       return assignment;
656:     };
659:     static Obj<ALE::Mesh> distributeMesh(const Obj<ALE::Mesh>& serialMesh, const int height = 0, const std::string& partitioner = "chaco") {
660:       MPI_Comm                          comm          = serialMesh->comm();
661:       const int                         dim           = serialMesh->getDimension();
662:       Obj<ALE::Mesh>                    parallelMesh  = new ALE::Mesh(comm, dim, serialMesh->debug());
663:       const Obj<ALE::Mesh::sieve_type>& parallelSieve = new ALE::Mesh::sieve_type(comm, serialMesh->debug());

665:       ALE_LOG_EVENT_BEGIN;
666:       parallelMesh->setSieve(parallelSieve);
667:       if (serialMesh->debug()) {serialMesh->view("Serial mesh");}

669:       // Distribute cones
670:       Obj<send_overlap_type> sendOverlap = new send_overlap_type(comm, serialMesh->debug());
671:       Obj<recv_overlap_type> recvOverlap = new recv_overlap_type(comm, serialMesh->debug());
672:       scatterBundle(serialMesh, dim, parallelMesh, sendOverlap, recvOverlap, partitioner, height);
673:       parallelMesh->setDistSendOverlap(sendOverlap);
674:       parallelMesh->setDistRecvOverlap(recvOverlap);

676:       // Distribute labels
677:       const typename bundle_type::labels_type& labels = serialMesh->getLabels();

679:       for(typename bundle_type::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
680:         if (parallelMesh->hasLabel(l_iter->first)) continue;
681:         const Obj<typename bundle_type::label_type>& serialLabel   = l_iter->second;
682:         const Obj<typename bundle_type::label_type>& parallelLabel = parallelMesh->createLabel(l_iter->first);
683:         // Create local label
684: #define NEW_LABEL
685: #ifdef NEW_LABEL
686:         parallelLabel->add(serialLabel, parallelSieve);
687: #else
688:         const Obj<typename bundle_type::label_type::traits::baseSequence>& base = serialLabel->base();

690:         for(typename bundle_type::label_type::traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
691:           if (parallelSieve->capContains(*b_iter) || parallelSieve->baseContains(*b_iter)) {
692:             parallelLabel->addArrow(*serialLabel->cone(*b_iter)->begin(), *b_iter);
693:           }
694:         }
695: #endif
696:         // Get remote labels
697:         sieveCompletion::scatterCones(serialLabel, parallelLabel, sendOverlap, recvOverlap);
698:       }

700:       // Distribute sections
701:       Obj<std::set<std::string> > sections = serialMesh->getRealSections();

703:       for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
704:         parallelMesh->setRealSection(*name, distributeSection(serialMesh->getRealSection(*name), parallelMesh, sendOverlap, recvOverlap));
705:       }
706:       sections = serialMesh->getIntSections();
707:       for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
708:         parallelMesh->setIntSection(*name, distributeSection(serialMesh->getIntSection(*name), parallelMesh, sendOverlap, recvOverlap));
709:       }
710:       sections = serialMesh->getArrowSections();

712:       for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
713:         parallelMesh->setArrowSection(*name, distributeArrowSection(serialMesh->getArrowSection(*name), serialMesh, parallelMesh, sendOverlap, recvOverlap));
714:       }
715:       if (parallelMesh->debug()) {parallelMesh->view("Parallel Mesh");}
716:       ALE_LOG_EVENT_END;
717:       return parallelMesh;
718:     };
721:     template<typename Section>
722:     static void updateSectionLocal(const Obj<Section>& oldSection, const Obj<bundle_type>& newBundle, const Obj<Section>& newSection) {
723:       const Obj<typename bundle_type::sieve_type>&    newSieve = newBundle->getSieve();
724:       const typename Section::atlas_type::chart_type& oldChart = oldSection->getChart();

726:       for(typename Section::atlas_type::chart_type::const_iterator c_iter = oldChart.begin(); c_iter != oldChart.end(); ++c_iter) {
727:         if (newSieve->capContains(*c_iter) || newSieve->baseContains(*c_iter)) {
728:           newSection->setFiberDimension(*c_iter, oldSection->getFiberDimension(*c_iter));
729:         }
730:       }
731:       newBundle->allocate(newSection);
732:       const typename Section::atlas_type::chart_type& newChart = newSection->getChart();

734:       for(typename Section::atlas_type::chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
735:         newSection->updatePointAll(*c_iter, oldSection->restrictPoint(*c_iter));
736:       }
737:     };
740:     template<typename RecvSection, typename Section>
741:     static void updateSectionRemote(const Obj<recv_overlap_type>& recvOverlap, const Obj<RecvSection>& recvSection, const Obj<bundle_type>& newBundle, const Obj<Section>& newSection) {
742:       Obj<typename recv_overlap_type::traits::baseSequence> recvPoints = recvOverlap->base();

744:       for(typename recv_overlap_type::traits::baseSequence::iterator r_iter = recvPoints->begin(); r_iter != recvPoints->end(); ++r_iter) {
745:         const Obj<typename recv_overlap_type::traits::coneSequence>&     recvPatches = recvOverlap->cone(*r_iter);
746:         const typename recv_overlap_type::traits::coneSequence::iterator end         = recvPatches->end();

748:         for(typename recv_overlap_type::traits::coneSequence::iterator p_iter = recvPatches->begin(); p_iter != end; ++p_iter) {
749:           newSection->addPoint(*r_iter, recvSection->getSection(*p_iter)->getFiberDimension(*r_iter));
750:         }
751:       }
752:       newBundle->reallocate(newSection);
753:       for(typename recv_overlap_type::traits::baseSequence::iterator r_iter = recvPoints->begin(); r_iter != recvPoints->end(); ++r_iter) {
754:         const Obj<typename recv_overlap_type::traits::coneSequence>&     recvPatches = recvOverlap->cone(*r_iter);
755:         const typename recv_overlap_type::traits::coneSequence::iterator end         = recvPatches->end();

757:         for(typename recv_overlap_type::traits::coneSequence::iterator p_iter = recvPatches->begin(); p_iter != end; ++p_iter) {
758:           if (recvSection->getSection(*p_iter)->getFiberDimension(*r_iter)) {
759:             newSection->updatePointAll(*r_iter, recvSection->getSection(*p_iter)->restrictPoint(*r_iter));
760:           }
761:         }
762:       }
763:     };
766:     template<typename Section>
767:     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) {
768:       if (serialSection->debug()) {
769:         serialSection->view("Serial Section");
770:       }
771:       typedef typename alloc_type::template rebind<typename Section::value_type>::other value_alloc_type;
772:       typedef ALE::Field<send_overlap_type, int, ALE::Section<point_type, typename Section::value_type, value_alloc_type> > send_section_type;
773:       typedef ALE::Field<recv_overlap_type, int, ALE::Section<point_type, typename Section::value_type, value_alloc_type> > recv_section_type;
774:       typedef ALE::New::SizeSection<Section> SectionSizer;
775:       Obj<Section>                 parallelSection = new Section(serialSection->comm(), serialSection->debug());
776:       const Obj<send_section_type> sendSection     = new send_section_type(serialSection->comm(), serialSection->debug());
777:       const Obj<recv_section_type> recvSection     = new recv_section_type(serialSection->comm(), sendSection->getTag(), serialSection->debug());
778:       const Obj<SectionSizer>      sizer           = new SectionSizer(serialSection);

780:       updateSectionLocal(serialSection, parallelBundle, parallelSection);
781:       sectionCompletion::completeSection(sendOverlap, recvOverlap, sizer, serialSection, sendSection, recvSection);
782:       updateSectionRemote(recvOverlap, recvSection, parallelBundle, parallelSection);
783:       if (parallelSection->debug()) {
784:         parallelSection->view("Parallel Section");
785:       }
786:       return parallelSection;
787:     };
790:     template<typename Section>
791:     static void updateArrowSectionLocal(const Obj<Section>& oldSection, const Obj<bundle_type>& newBundle, const Obj<Section>& newSection) {
792:       const Obj<typename bundle_type::sieve_type>&    newSieve = newBundle->getSieve();
793:       const typename Section::atlas_type::chart_type& oldChart = oldSection->getChart();

795:       for(typename Section::atlas_type::chart_type::const_iterator c_iter = oldChart.begin(); c_iter != oldChart.end(); ++c_iter) {
796:         // Dmitry should provide a Sieve::contains(MinimalArrow) method
797:         if (newSieve->capContains(c_iter->source) && newSieve->baseContains(c_iter->target)) {
798:           newSection->setFiberDimension(*c_iter, oldSection->getFiberDimension(*c_iter));
799:         }
800:       }
801:       //newBundle->allocate(newSection);
802:       const typename Section::atlas_type::chart_type& newChart = newSection->getChart();

804:       for(typename Section::atlas_type::chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
805:         newSection->updatePointAll(*c_iter, oldSection->restrictPoint(*c_iter));
806:       }
807:     };
810:     template<typename RecvSection, typename Section>
811:     static void updateArrowSectionRemote(const Obj<recv_overlap_type>& recvOverlap, const Obj<RecvSection>& recvSection, const Obj<bundle_type>& newBundle, const Obj<Section>& newSection) {
812:       Obj<typename recv_overlap_type::traits::baseSequence> recvPoints = recvOverlap->base();

814:       for(typename recv_overlap_type::traits::baseSequence::iterator r_iter = recvPoints->begin(); r_iter != recvPoints->end(); ++r_iter) {
815:         const Obj<typename bundle_type::sieve_type::traits::coneSequence>&     cone = newBundle->getSieve()->cone(*r_iter);
816:         const typename bundle_type::sieve_type::traits::coneSequence::iterator end  = cone->end();

818:         for(typename bundle_type::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
819:           newSection->setFiberDimension(typename Section::point_type(*c_iter, *r_iter), 1);
820:         }
821:       }
822:       //newBundle->reallocate(newSection);
823:       for(typename recv_overlap_type::traits::baseSequence::iterator r_iter = recvPoints->begin(); r_iter != recvPoints->end(); ++r_iter) {
824:         const Obj<typename recv_overlap_type::traits::coneSequence>&     recvPatches = recvOverlap->cone(*r_iter);
825:         const typename recv_overlap_type::traits::coneSequence::iterator recvEnd     = recvPatches->end();

827:         for(typename recv_overlap_type::traits::coneSequence::iterator p_iter = recvPatches->begin(); p_iter != recvEnd; ++p_iter) {
828:           const Obj<typename RecvSection::section_type>& section = recvSection->getSection(*p_iter);

830:           if (section->getFiberDimension(*r_iter)) {
831:             const Obj<typename bundle_type::sieve_type::traits::coneSequence>&     cone    = newBundle->getSieve()->cone(*r_iter);
832:             const typename bundle_type::sieve_type::traits::coneSequence::iterator end     = cone->end();
833:             const typename RecvSection::value_type                                *values  = section->restrictPoint(*r_iter);
834:             int                                                                    c       = -1;

836:             for(typename bundle_type::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
837:               newSection->updatePoint(typename Section::point_type(*c_iter, *r_iter), &values[++c]);
838:             }
839:           }
840:         }
841:       }
842:     };
845:     template<typename Section>
846:     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) {
847:       if (serialSection->debug()) {
848:         serialSection->view("Serial ArrowSection");
849:       }
850:       typedef typename alloc_type::template rebind<typename Section::value_type>::other value_alloc_type;
851:       typedef ALE::Field<send_overlap_type, int, ALE::Section<point_type, typename Section::value_type, value_alloc_type> > send_section_type;
852:       typedef ALE::Field<recv_overlap_type, int, ALE::Section<point_type, typename Section::value_type, value_alloc_type> > recv_section_type;
853:       typedef ALE::New::ConeSizeSection<bundle_type, sieve_type> SectionSizer;
854:       typedef ALE::New::ArrowSection<sieve_type, Section>        ArrowFiller;
855:       Obj<Section>                 parallelSection = new Section(serialSection->comm(), serialSection->debug());
856:       const Obj<send_section_type> sendSection     = new send_section_type(serialSection->comm(), serialSection->debug());
857:       const Obj<recv_section_type> recvSection     = new recv_section_type(serialSection->comm(), sendSection->getTag(), serialSection->debug());
858:       const Obj<SectionSizer>      sizer           = new SectionSizer(serialBundle, serialBundle->getSieve());
859:       const Obj<ArrowFiller>       filler          = new ArrowFiller(serialBundle->getSieve(), serialSection);

861:       updateArrowSectionLocal(serialSection, parallelBundle, parallelSection);
862:       sectionCompletion::completeSection(sendOverlap, recvOverlap, sizer, filler, sendSection, recvSection);
863:       updateArrowSectionRemote(recvOverlap, recvSection, parallelBundle, parallelSection);
864:       if (parallelSection->debug()) {
865:         parallelSection->view("Parallel ArrowSection");
866:       }
867:       return parallelSection;
868:     };
869:     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) {
870:       typedef int part_type;
871:       const Obj<sieve_type>& sieve    = bundle->getSieve();
872:       const Obj<sieve_type>& sieveNew = bundleNew->getSieve();
873:       const int              rank     = bundle->commRank();
874:       const int              debug    = bundle->debug();

876:       // 1) Form partition point overlap a priori
877:       if (rank == 0) {
878:         for(int p = 1; p < sieve->commSize(); p++) {
879:           // The arrow is from remote partition point 0 on rank p to local partition point 0
880:           recvOverlap->addCone(p, 0, 0);
881:         }
882:       } else {
883:         // The arrow is from local partition point 0 to remote partition point 0 on rank 0
884:         sendOverlap->addCone(0, 0, 0);
885:       }
886:       if (debug) {
887:         sendOverlap->view("Send overlap for partition");
888:         recvOverlap->view("Receive overlap for partition");
889:       }
890:       // 2) Partition the mesh
891:       int        numCells   = bundle->heightStratum(0)->size();
892:       part_type *assignment = new part_type[numCells];

894:       for(int c = 0; c < numCells; ++c) {
895:         assignment[c] = 0;
896:       }
897:       // 3) Scatter the sieve
898:       sieveCompletion::scatterSieve(bundle, sieve, dim, sieveNew, sendOverlap, recvOverlap, 0, numCells, assignment);
899:       bundleNew->stratify();
900:       // 4) Cleanup
901:       if (assignment != NULL) delete [] assignment;
902:     };
905:     static Obj<ALE::Mesh> unifyMesh(const Obj<ALE::Mesh>& parallelMesh) {
906:       const int                         dim         = parallelMesh->getDimension();
907:       Obj<ALE::Mesh>                    serialMesh  = new ALE::Mesh(parallelMesh->comm(), dim, parallelMesh->debug());
908:       const Obj<ALE::Mesh::sieve_type>& serialSieve = new ALE::Mesh::sieve_type(parallelMesh->comm(), parallelMesh->debug());

910:       ALE_LOG_EVENT_BEGIN;
911:       serialMesh->setSieve(serialSieve);
912:       if (parallelMesh->debug()) {
913:         parallelMesh->view("Parallel topology");
914:       }

916:       // Unify cones
917:       Obj<send_overlap_type> sendOverlap = new send_overlap_type(serialMesh->comm(), serialMesh->debug());
918:       Obj<recv_overlap_type> recvOverlap = new recv_overlap_type(serialMesh->comm(), serialMesh->debug());
919:       unifyBundle(parallelMesh, dim, serialMesh, sendOverlap, recvOverlap);
920:       serialMesh->setDistSendOverlap(sendOverlap);
921:       serialMesh->setDistRecvOverlap(recvOverlap);

923:       // Unify labels
924:       const typename bundle_type::labels_type& labels = parallelMesh->getLabels();

926:       for(typename bundle_type::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
927:         if (serialMesh->hasLabel(l_iter->first)) continue;
928:         const Obj<typename bundle_type::label_type>& parallelLabel = l_iter->second;
929:         const Obj<typename bundle_type::label_type>& serialLabel   = serialMesh->createLabel(l_iter->first);

931:         sieveCompletion::scatterCones(parallelLabel, serialLabel, sendOverlap, recvOverlap);
932:       }

934:       // Unify coordinates
935:       Obj<std::set<std::string> > sections = parallelMesh->getRealSections();

937:       for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
938:         serialMesh->setRealSection(*name, distributeSection(parallelMesh->getRealSection(*name), serialMesh, sendOverlap, recvOverlap));
939:       }
940:       sections = parallelMesh->getIntSections();
941:       for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
942:         serialMesh->setIntSection(*name, distributeSection(parallelMesh->getIntSection(*name), serialMesh, sendOverlap, recvOverlap));
943:       }
944:       sections = parallelMesh->getArrowSections();
945:       for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
946:         serialMesh->setArrowSection(*name, distributeArrowSection(parallelMesh->getArrowSection(*name), parallelMesh, serialMesh, sendOverlap, recvOverlap));
947:       }
948:       if (serialMesh->debug()) {serialMesh->view("Serial Mesh");}
949:       ALE_LOG_EVENT_END;
950:       return serialMesh;
951:     };
952:   public: // Do not like these
955:     // This is just crappy. We could introduce another phase to find out exactly what
956:     //   indices people do not have in the global order after communication
957:     template<typename OrigSendOverlap, typename OrigRecvOverlap, typename SendSection, typename RecvSection>
958:     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) {
959:       const typename SendSection::sheaf_type& sendRanks = sendSection->getPatches();
960:       const typename RecvSection::sheaf_type& recvRanks = recvSection->getPatches();

962:       for(typename SendSection::sheaf_type::const_iterator p_iter = sendRanks.begin(); p_iter != sendRanks.end(); ++p_iter) {
963:         const typename SendSection::patch_type&               rank    = p_iter->first;
964:         const Obj<typename SendSection::section_type>&        section = p_iter->second;
965:         const typename SendSection::section_type::chart_type& chart   = section->getChart();

967:         for(typename SendSection::section_type::chart_type::const_iterator b_iter = chart.begin(); b_iter != chart.end(); ++b_iter) {
968:           const typename SendSection::value_type *points = section->restrictPoint(*b_iter);
969:           const int                               size   = section->getFiberDimension(*b_iter);

971:           for(int p = 0; p < size; p++) {
972:             if (origSendOverlap->support(points[p])->size() == 0) {
973:               sendOverlap->addArrow(points[p], rank, points[p]);
974:             }
975:           }
976:         }
977:       }
978:       for(typename RecvSection::sheaf_type::const_iterator p_iter = recvRanks.begin(); p_iter != recvRanks.end(); ++p_iter) {
979:         const typename RecvSection::patch_type&               rank    = p_iter->first;
980:         const Obj<typename RecvSection::section_type>&        section = p_iter->second;
981:         const typename RecvSection::section_type::chart_type& chart   = section->getChart();

983:         for(typename RecvSection::section_type::chart_type::const_iterator b_iter = chart.begin(); b_iter != chart.end(); ++b_iter) {
984:           const typename RecvSection::value_type *points = section->restrictPoint(*b_iter);
985:           const int                               size   = section->getFiberDimension(*b_iter);

987:           for(int p = 0; p < size; p++) {
988:             if (origRecvOverlap->support(rank, points[p])->size() == 0) {
989:               recvOverlap->addArrow(rank, points[p], points[p]);
990:             }
991:           }
992:         }
993:       }
994:     };
997:     template<typename RecvOverlap, typename RecvSection>
998:     static void updateSieve(const Obj<RecvOverlap>& recvOverlap, const Obj<RecvSection>& recvSection, const Obj<sieve_type>& sieve) {
999: #if 1
1000:       Obj<typename RecvOverlap::traits::baseSequence> recvPoints = recvOverlap->base();

1002:       for(typename RecvOverlap::traits::baseSequence::iterator p_iter = recvPoints->begin(); p_iter != recvPoints->end(); ++p_iter) {
1003:         const Obj<typename RecvOverlap::traits::coneSequence>& ranks      = recvOverlap->cone(*p_iter);
1004:         const typename RecvOverlap::target_type&               localPoint = *p_iter;

1006:         for(typename RecvOverlap::traits::coneSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
1007:           const typename RecvOverlap::target_type&       remotePoint = r_iter.color();
1008:           const int                                      rank        = *r_iter;
1009:           const Obj<typename RecvSection::section_type>& section     = recvSection->getSection(rank);
1010:           const typename RecvSection::value_type        *points      = section->restrictPoint(remotePoint);
1011:           const int                                      size        = section->getFiberDimension(remotePoint);
1012:           int                                            c           = 0;

1014:           ///std::cout << "["<<recvSection->commRank()<<"]: Receiving " << size << " points from rank " << rank << std::endl;
1015:           for(int p = 0; p < size; p++) {
1016:             // rank -- remote point --> local point
1017:             if (recvOverlap->support(rank, points[p])->size()) {
1018:               sieve->addArrow(*recvOverlap->support(rank, points[p])->begin(), localPoint, c);
1019:               ///std::cout << "["<<recvSection->commRank()<<"]:   1Adding arrow " << *recvOverlap->support(rank, points[p])->begin() << "("<<points[p]<<") --> " << localPoint << std::endl;
1020:             } else {
1021:               sieve->addArrow(points[p], localPoint, c);
1022:               ///std::cout << "["<<recvSection->commRank()<<"]:   2Adding arrow " << points[p] << " --> " << localPoint << std::endl;
1023:             }
1024:           }
1025:         }
1026:       }
1027: #else
1028:       const typename RecvSection::sheaf_type& ranks = recvSection->getPatches();

1030:       for(typename RecvSection::sheaf_type::const_iterator p_iter = ranks.begin(); p_iter != ranks.end(); ++p_iter) {
1031:         const Obj<typename RecvSection::section_type>&        section = p_iter->second;
1032:         const typename RecvSection::section_type::chart_type& chart   = section->getChart();

1034:         for(typename RecvSection::section_type::chart_type::const_iterator b_iter = chart.begin(); b_iter != chart.end(); ++b_iter) {
1035:           const typename RecvSection::value_type *points = section->restrictPoint(*b_iter);
1036:           int                                     size   = section->getFiberDimension(*b_iter);
1037:           int                                     c      = 0;

1039:           std::cout << "["<<recvSection->commRank()<<"]: Receiving " << size << " points from rank " << p_iter->first << std::endl;
1040:           for(int p = 0; p < size; p++) {
1041:             //sieve->addArrow(points[p], *b_iter, c++);
1042:             sieve->addArrow(points[p], *b_iter, c);
1043:             std::cout << "["<<recvSection->commRank()<<"]:   Adding arrow " << points[p] << " --> " << *b_iter << std::endl;
1044:           }
1045:         }
1046:       }
1047: #endif
1048:     };
1051:     template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
1052:     static void coneCompletion(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<bundle_type>& bundle, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
1053:       if (sendOverlap->commSize() == 1) return;
1054:       // Distribute cones
1055:       const Obj<sieve_type>&                                 sieve           = bundle->getSieve();
1056:       const Obj<typename sieveCompletion::topology_type>     secTopology     = sieveCompletion::completion::createSendTopology(sendOverlap);
1057:       const Obj<typename sieveCompletion::cone_size_section> coneSizeSection = new typename sieveCompletion::cone_size_section(bundle, sieve);
1058:       const Obj<typename sieveCompletion::cone_section>      coneSection     = new typename sieveCompletion::cone_section(sieve);
1059:       sieveCompletion::completion::completeSection(sendOverlap, recvOverlap, coneSizeSection, coneSection, sendSection, recvSection);
1060:       // Update cones
1061:       updateSieve(recvOverlap, recvSection, sieve);
1062:     };
1065:     template<typename Section>
1066:     static void completeSection(const Obj<bundle_type>& bundle, const Obj<Section>& section) {
1067:       typedef typename Distribution<bundle_type>::sieveCompletion sieveCompletion;
1068:       typedef typename bundle_type::send_overlap_type             send_overlap_type;
1069:       typedef typename bundle_type::recv_overlap_type             recv_overlap_type;
1070:       typedef typename Section::value_type                        value_type;
1071:       typedef typename alloc_type::template rebind<typename Section::value_type>::other value_alloc_type;
1072:       typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
1073:       typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
1074:       typedef ALE::New::SizeSection<Section>                                SectionSizer;
1075:       const int debug = section->debug();

1077:       bundle->constructOverlap();
1078:       const Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
1079:       const Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();
1080:       const Obj<send_section_type> sendSection = new send_section_type(section->comm(), section->debug());
1081:       const Obj<recv_section_type> recvSection = new recv_section_type(section->comm(), sendSection->getTag(), section->debug());
1082:       const Obj<SectionSizer>      sizer       = new SectionSizer(section);

1084:       sectionCompletion::completeSection(sendOverlap, recvOverlap, sizer, section, sendSection, recvSection);
1085:       // Update section with remote data
1086:       const Obj<typename recv_overlap_type::traits::baseSequence> recvPoints = bundle->getRecvOverlap()->base();

1088:       for(typename recv_overlap_type::traits::baseSequence::iterator r_iter = recvPoints->begin(); r_iter != recvPoints->end(); ++r_iter) {
1089:         const Obj<typename recv_overlap_type::traits::coneSequence>&     recvPatches = recvOverlap->cone(*r_iter);
1090:         const typename recv_overlap_type::traits::coneSequence::iterator end         = recvPatches->end();

1092:         for(typename recv_overlap_type::traits::coneSequence::iterator p_iter = recvPatches->begin(); p_iter != end; ++p_iter) {
1093:           if (recvSection->getSection(*p_iter)->getFiberDimension(p_iter.color())) {
1094:             if (debug) {std::cout << "["<<section->commRank()<<"]Completed point " << *r_iter << std::endl;}
1095:             section->updateAddPoint(*r_iter, recvSection->getSection(*p_iter)->restrictPoint(p_iter.color()));
1096:           }
1097:         }
1098:       }
1099:     };
1100:   };
1101: }
1102: #endif