Actual source code: Distribution.hh

  1: #ifndef included_ALE_Distribution_hh
  2: #define included_ALE_Distribution_hh

  4: #ifndef  included_ALE_Mesh_hh
  5: #include <Mesh.hh>
  6: #endif

  8: #ifndef  included_ALE_Completion_hh
  9: #include <Completion.hh>
 10: #endif

 12: // Attempt to unify all of the distribution mechanisms:
 13: //   one to many  (distributeMesh)
 14: //   many to one  (unifyMesh)
 15: //   many to many (Numbering)
 16: // as well as things being distributed
 17: //   Section
 18: //   Sieve        (This sends two sections, the points and cones)
 19: //   Numbering    (Should be an integer section)
 20: //   Global Order (should be an integer section with extra methods)
 21: //
 22: // 0) Create the new object to hold the communicated data
 23: //
 24: // 1) Create Overlap
 25: //    There may be special ways to do this based upon what we know at the time
 26: //
 27: // 2) Create send and receive sections over the interface
 28: //    These have a flat topology now, consisting only of the overlap nodes
 29: //    We could make a full topology on the overlap (maybe it is necessary for higher order)
 30: //
 31: // 3) Communication section
 32: //    Create sizer sections on interface (uses constant sizer)
 33: //    Communicate sizes on interface (uses custom filler)
 34: //      Fill send section
 35: //      sendSection->startCommunication();
 36: //      recvSection->startCommunication();
 37: //      sendSection->endCommunication();
 38: //      recvSection->endCommunication();
 39: //
 40: //    Create section on interface (uses previous sizer)
 41: //    Communicate values on interface (uses custom filler)
 42: //      Same stuff as above
 43: //
 44: // 4) Update new section with old local values (can be done in between the communication?)
 45: //    Loop over patches in new topology
 46: //      Loop over chart from patch in old atlas
 47: //        If this point is in the new sieve from patch
 48: //          Set to old fiber dimension
 49: //    Order and allocate new section
 50: //    Repeat loop, but update values
 51: //
 52: // 5) Update new section with old received values
 53: //    Loop over patches in discrete topology of receive section (these are ranks)
 54: //      Loop over base of discrete sieve (we should transform this to a chart to match above)
 55: //        Get new patch from overlap, or should the receive patches be <rank, patch>?
 56: //        Guaranteed to be in the new sieve from patch (but we could check anyway)
 57: //          Set to recevied fiber dimension
 58: //    Order and allocate new section
 59: //    Repeat loop, but update values
 60: //
 61: // 6) Synchronize PETSc tags (can I get around this?)
 62: namespace ALE {
 63:   template<typename Mesh, typename Partitioner = ALE::Partitioner<> >
 64:   class DistributionNew {
 65:   public:
 66:     typedef Partitioner                                   partitioner_type;
 67:     typedef typename Mesh::point_type                     point_type;
 68:     typedef OrientedPoint<point_type>                     oriented_point_type;
 69:     typedef typename Partitioner::part_type               rank_type;
 70:     typedef ALE::ISection<rank_type, point_type>          partition_type;
 71:     typedef ALE::Section<ALE::Pair<int, point_type>, point_type>          cones_type;
 72:     typedef ALE::Section<ALE::Pair<int, point_type>, oriented_point_type> oriented_cones_type;
 73:   public:
 74:     template<typename Sieve, typename NewSieve, typename Renumbering, typename SendOverlap, typename RecvOverlap>
 75:     static Obj<cones_type> completeCones(const Obj<Sieve>& sieve, const Obj<NewSieve>& newSieve, Renumbering& renumbering, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap) {
 76:       typedef ALE::ConeSection<Sieve> cones_wrapper_type;
 77:       Obj<cones_wrapper_type> cones        = new cones_wrapper_type(sieve);
 78:       Obj<cones_type>         overlapCones = new cones_type(sieve->comm(), sieve->debug());

 80:       ALE::Pullback::SimpleCopy::copy(sendMeshOverlap, recvMeshOverlap, cones, overlapCones);
 81:       if (sieve->debug()) {overlapCones->view("Overlap Cones");}
 82:       // Inserts cones into parallelMesh (must renumber here)
 83:       ALE::Pullback::InsertionBinaryFusion::fuse(overlapCones, recvMeshOverlap, renumbering, newSieve);
 84:       return overlapCones;
 85:     };
 86:     template<typename Sieve, typename NewSieve, typename SendOverlap, typename RecvOverlap>
 87:     static Obj<oriented_cones_type> completeConesV(const Obj<Sieve>& sieve, const Obj<NewSieve>& newSieve, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap) {
 88:       typedef ALE::OrientedConeSectionV<Sieve> oriented_cones_wrapper_type;
 89:       Obj<oriented_cones_wrapper_type> cones        = new oriented_cones_wrapper_type(sieve);
 90:       Obj<oriented_cones_type>         overlapCones = new oriented_cones_type(sieve->comm(), sieve->debug());

 92:       ALE::Pullback::SimpleCopy::copy(sendMeshOverlap, recvMeshOverlap, cones, overlapCones);
 93:       if (sieve->debug()) {overlapCones->view("Overlap Cones");}
 94:       ALE::Pullback::InsertionBinaryFusion::fuse(overlapCones, recvMeshOverlap, newSieve);
 95:       return overlapCones;
 96:     };
 97:     template<typename Sieve, typename NewSieve, typename Renumbering, typename SendOverlap, typename RecvOverlap>
 98:     static Obj<oriented_cones_type> completeConesV(const Obj<Sieve>& sieve, const Obj<NewSieve>& newSieve, Renumbering& renumbering, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap) {
 99:       typedef ALE::OrientedConeSectionV<Sieve> oriented_cones_wrapper_type;
100:       Obj<oriented_cones_wrapper_type> cones        = new oriented_cones_wrapper_type(sieve);
101:       Obj<oriented_cones_type>         overlapCones = new oriented_cones_type(sieve->comm(), sieve->debug());

103:       ALE::Pullback::SimpleCopy::copy(sendMeshOverlap, recvMeshOverlap, cones, overlapCones);
104:       if (sieve->debug()) {overlapCones->view("Overlap Cones");}
105:       // Inserts cones into parallelMesh (must renumber here)
106:       ALE::Pullback::InsertionBinaryFusion::fuse(overlapCones, recvMeshOverlap, renumbering, newSieve);
107:       return overlapCones;
108:     };
109:     // Given a partition of sieve points, copy the mesh pieces to each process and fuse into the new mesh
110:     //   Return overlaps for the cone communication
111:     template<typename Renumbering, typename NewMesh, typename SendOverlap, typename RecvOverlap>
112:     static void completeMesh(const Obj<Mesh>& mesh, const Obj<partition_type>& partition, Renumbering& renumbering, const Obj<NewMesh>& newMesh, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap) {
113:       typedef ALE::Sifter<rank_type,rank_type,rank_type> part_send_overlap_type;
114:       typedef ALE::Sifter<rank_type,rank_type,rank_type> part_recv_overlap_type;
115:       const Obj<part_send_overlap_type> sendOverlap = new part_send_overlap_type(partition->comm());
116:       const Obj<part_recv_overlap_type> recvOverlap = new part_recv_overlap_type(partition->comm());

118:       // Create overlap for partition points
119:       //   TODO: This needs to be generalized for multiple sources
120:       Partitioner::createDistributionPartOverlap(sendOverlap, recvOverlap);
121:       // Communicate partition pieces to processes
122:       Obj<partition_type> overlapPartition = new partition_type(partition->comm(), partition->debug());

124:       overlapPartition->setChart(partition->getChart());
125:       ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, partition, overlapPartition);
126:       // Create renumbering
127:       const int         rank           = partition->commRank();
128:       const point_type *localPoints    = partition->restrictPoint(rank);
129:       const int         numLocalPoints = partition->getFiberDimension(rank);

131:       for(point_type p = 0; p < numLocalPoints; ++p) {
132:         renumbering[localPoints[p]] = p;
133:       }
134:       const Obj<typename part_recv_overlap_type::traits::baseSequence> rPoints    = recvOverlap->base();
135:       point_type                                                       localPoint = numLocalPoints;

137:       for(typename part_recv_overlap_type::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rPoints->end(); ++p_iter) {
138:         const Obj<typename part_recv_overlap_type::coneSequence>& ranks           = recvOverlap->cone(*p_iter);
139:         const rank_type&                                          remotePartPoint = ranks->begin().color();
140:         const point_type                                         *points          = overlapPartition->restrictPoint(remotePartPoint);
141:         const int                                                 numPoints       = overlapPartition->getFiberDimension(remotePartPoint);

143:         for(int i = 0; i < numPoints; ++i) {
144:           renumbering[points[i]] = localPoint++;
145:         }
146:       }
147:       // Create mesh overlap from partition overlap
148:       //   TODO: Generalize to redistribution (receive from multiple sources)
149:       Partitioner::createDistributionMeshOverlap(partition, recvOverlap, renumbering, overlapPartition, sendMeshOverlap, recvMeshOverlap);
150:       // Send cones
151:       completeCones(mesh->getSieve(), newMesh->getSieve(), renumbering, sendMeshOverlap, recvMeshOverlap);
152:     };
153:     template<typename Renumbering, typename NewMesh, typename SendOverlap, typename RecvOverlap>
154:     static void completeBaseV(const Obj<Mesh>& mesh, const Obj<partition_type>& partition, Renumbering& renumbering, const Obj<NewMesh>& newMesh, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap) {
155:       typedef ALE::Sifter<rank_type,rank_type,rank_type> part_send_overlap_type;
156:       typedef ALE::Sifter<rank_type,rank_type,rank_type> part_recv_overlap_type;
157:       const Obj<part_send_overlap_type> sendOverlap = new part_send_overlap_type(partition->comm());
158:       const Obj<part_recv_overlap_type> recvOverlap = new part_recv_overlap_type(partition->comm());

160:       // Create overlap for partition points
161:       //   TODO: This needs to be generalized for multiple sources
162:       Partitioner::createDistributionPartOverlap(sendOverlap, recvOverlap);
163:       // Communicate partition pieces to processes
164:       Obj<partition_type> overlapPartition = new partition_type(partition->comm(), partition->debug());

166:       overlapPartition->setChart(partition->getChart());
167:       ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, partition, overlapPartition);
168:       // Create renumbering
169:       const int         rank           = partition->commRank();
170:       const point_type *localPoints    = partition->restrictPoint(rank);
171:       const int         numLocalPoints = partition->getFiberDimension(rank);

173:       for(point_type p = 0; p < numLocalPoints; ++p) {
174:         ///std::cout <<"["<<partition->commRank()<<"]: local renumbering " << localPoints[p] << " --> " << p << std::endl;
175:         renumbering[localPoints[p]] = p;
176:       }
177:       const Obj<typename part_recv_overlap_type::traits::baseSequence> rPoints    = recvOverlap->base();
178:       point_type                                                       localPoint = numLocalPoints;

180:       for(typename part_recv_overlap_type::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rPoints->end(); ++p_iter) {
181:         const Obj<typename part_recv_overlap_type::coneSequence>& ranks           = recvOverlap->cone(*p_iter);
182:         const rank_type&                                          remotePartPoint = ranks->begin().color();
183:         const point_type                                         *points          = overlapPartition->restrictPoint(remotePartPoint);
184:         const int                                                 numPoints       = overlapPartition->getFiberDimension(remotePartPoint);

186:         for(int i = 0; i < numPoints; ++i) {
187:           ///std::cout <<"["<<partition->commRank()<<"]: remote renumbering " << points[i] << " --> " << localPoint << std::endl;
188:           renumbering[points[i]] = localPoint++;
189:         }
190:       }
191:       newMesh->getSieve()->setChart(typename NewMesh::sieve_type::chart_type(0, renumbering.size()));
192:       // Create mesh overlap from partition overlap
193:       //   TODO: Generalize to redistribution (receive from multiple sources)
194:       Partitioner::createDistributionMeshOverlap(partition, recvOverlap, renumbering, overlapPartition, sendMeshOverlap, recvMeshOverlap);
195:     };
196:     template<typename NewMesh, typename Renumbering, typename SendOverlap, typename RecvOverlap>
197:     static Obj<partition_type> distributeMesh(const Obj<Mesh>& mesh, const Obj<NewMesh>& newMesh, Renumbering& renumbering, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap, const int height = 0) {
198:       const Obj<partition_type> cellPartition = new partition_type(mesh->comm(), 0, mesh->commSize(), mesh->debug());
199:       const Obj<partition_type> partition     = new partition_type(mesh->comm(), 0, mesh->commSize(), mesh->debug());

201:       // Create the cell partition
202:       Partitioner::createPartition(mesh, cellPartition, height);
203:       if (mesh->debug()) {
204:         PetscViewer    viewer;

207:         cellPartition->view("Cell Partition");
208:         PetscViewerCreate(mesh->comm(), &viewer);CHKERRXX(ierr);
209:         PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);CHKERRXX(ierr);
210:         PetscViewerFileSetName(viewer, "mesh.vtk");CHKERRXX(ierr);
211:         ///TODO MeshView_Sieve_Ascii(mesh, cellPartition, viewer);CHKERRXX(ierr);
212:         PetscViewerDestroy(viewer);CHKERRXX(ierr);
213:       }
214:       // Close the partition over sieve points
215:       Partitioner::createPartitionClosure(mesh, cellPartition, partition, height);
216:       if (mesh->debug()) {partition->view("Partition");}
217:       // Create the remote meshes
218:       completeMesh(mesh, partition, renumbering, newMesh, sendMeshOverlap, recvMeshOverlap);
219:       // Create the local mesh
220:       Partitioner::createLocalMesh(mesh, partition, renumbering, newMesh, height);
221:       newMesh->stratify();
222:       return partition;
223:     };
224:     template<typename NewMesh, typename Renumbering, typename SendOverlap, typename RecvOverlap>
225:     static Obj<partition_type> distributeMeshAndSections(const Obj<Mesh>& mesh, const Obj<NewMesh>& newMesh, Renumbering& renumbering, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap, const int height = 0) {
226:       Obj<partition_type> partition = distributeMesh(mesh, newMesh, renumbering, sendMeshOverlap, recvMeshOverlap, height);

228:       // Distribute the coordinates
229:       const Obj<typename Mesh::real_section_type>& coordinates         = mesh->getRealSection("coordinates");
230:       const Obj<typename Mesh::real_section_type>& parallelCoordinates = newMesh->getRealSection("coordinates");

232:       newMesh->setupCoordinates(parallelCoordinates);
233:       distributeSection(coordinates, partition, renumbering, sendMeshOverlap, recvMeshOverlap, parallelCoordinates);
234:       // Distribute other sections
235:       if (mesh->getRealSections()->size() > 1) {
236:         Obj<std::set<std::string> > names = mesh->getRealSections();

238:         for(std::set<std::string>::const_iterator n_iter = names->begin(); n_iter != names->end(); ++n_iter) {
239:           if (*n_iter == "coordinates")   continue;
240:           distributeSection(mesh->getRealSection(*n_iter), partition, renumbering, sendMeshOverlap, recvMeshOverlap, newMesh->getRealSection(*n_iter));
241:         }
242:       }
243:       if (mesh->getIntSections()->size() > 0) {
244:         Obj<std::set<std::string> > names = mesh->getIntSections();

246:         for(std::set<std::string>::const_iterator n_iter = names->begin(); n_iter != names->end(); ++n_iter) {
247:           distributeSection(mesh->getIntSection(*n_iter), partition, renumbering, sendMeshOverlap, recvMeshOverlap, newMesh->getIntSection(*n_iter));
248:         }
249:       }
250:       if (mesh->getArrowSections()->size() > 1) {
251:         throw ALE::Exception("Need to distribute more arrow sections");
252:       }
253:       // Distribute labels
254:       const typename Mesh::labels_type& labels = mesh->getLabels();

256:       for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
257:         if (newMesh->hasLabel(l_iter->first)) continue;
258:         const Obj<typename Mesh::label_type>& origLabel = l_iter->second;
259:         const Obj<typename Mesh::label_type>& newLabel  = newMesh->createLabel(l_iter->first);
260:         // Get remote labels
261:         ALE::New::Completion<Mesh,typename Mesh::point_type>::scatterCones(origLabel, newLabel, sendMeshOverlap, recvMeshOverlap, renumbering);
262:         // Create local label
263:         newLabel->add(origLabel, newMesh->getSieve(), renumbering);
264:       }
265:       return partition;
266:     };
267:     template<typename NewMesh, typename Renumbering, typename SendOverlap, typename RecvOverlap>
268:     static Obj<partition_type> distributeMeshV(const Obj<Mesh>& mesh, const Obj<NewMesh>& newMesh, Renumbering& renumbering, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap, const int height = 0) {
269:       const Obj<partition_type> cellPartition = new partition_type(mesh->comm(), 0, mesh->commSize(), mesh->debug());
270:       const Obj<partition_type> partition     = new partition_type(mesh->comm(), 0, mesh->commSize(), mesh->debug());

272:       PETSc::Log::Event("DistributeMesh").begin();
273:       // Create the cell partition
274:       Partitioner::createPartitionV(mesh, cellPartition, height);
275:       if (mesh->debug()) {
276:         PetscViewer    viewer;

279:         cellPartition->view("Cell Partition");
280:         PetscViewerCreate(mesh->comm(), &viewer);CHKERRXX(ierr);
281:         PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);CHKERRXX(ierr);
282:         PetscViewerFileSetName(viewer, "mesh.vtk");CHKERRXX(ierr);
283:         ///TODO MeshView_Sieve_Ascii(mesh, cellPartition, viewer);CHKERRXX(ierr);
284:         PetscViewerDestroy(viewer);CHKERRXX(ierr);
285:       }
286:       // Close the partition over sieve points
287:       Partitioner::createPartitionClosureV(mesh, cellPartition, partition, height);
288:       if (mesh->debug()) {partition->view("Partition");}
289:       // Create the remote bases
290:       completeBaseV(mesh, partition, renumbering, newMesh, sendMeshOverlap, recvMeshOverlap);
291:       // Size the local mesh
292:       Partitioner::sizeLocalMeshV(mesh, partition, renumbering, newMesh, height);
293:       // Create the remote meshes
294:       completeConesV(mesh->getSieve(), newMesh->getSieve(), renumbering, sendMeshOverlap, recvMeshOverlap);
295:       // Create the local mesh
296:       Partitioner::createLocalMeshV(mesh, partition, renumbering, newMesh, height);
297:       newMesh->getSieve()->symmetrize();
298:       newMesh->stratify();
299:       PETSc::Log::Event("DistributeMesh").end();
300:       return partition;
301:     };
302:     // distributeMeshV:
303:     //   createPartitionV (can be dumb)
304:     //   createPartitionClosureV (should be low memory)
305:     //   completeBaseV (???)
306:     //     Partitioner::createDistributionPartOverlap (low memory)
307:     //     copy points to partitions (uses small overlap and fake sections)
308:     //     renumber (map is potentially big, can measure)
309:     //     Partitioner::createDistributionMeshOverlap (should be large for distribution)
310:     //       sendMeshOverlap is localPoint--- remotePoint --->remoteRank
311:     //       recvMeshOverlap is remoteRank--- remotePoint --->localPoint
312:     //   sizeLocalMeshV (should be low memory)
313:     //   completeConesV (???)
314:     //   createLocalMesh (should be low memory)
315:     //   symmetrize
316:     //   stratify
317:     template<typename NewMesh>
318:     static void distributeMeshAndSectionsV(const Obj<Mesh>& mesh, const Obj<NewMesh>& newMesh) {
319:       typedef typename Mesh::point_type point_type;

321:       const Obj<typename Mesh::send_overlap_type> sendMeshOverlap = new typename Mesh::send_overlap_type(mesh->comm(), mesh->debug());
322:       const Obj<typename Mesh::recv_overlap_type> recvMeshOverlap = new typename Mesh::recv_overlap_type(mesh->comm(), mesh->debug());
323:       std::map<point_type,point_type>&            renumbering     = newMesh->getRenumbering();
324:       // Distribute the mesh
325:       Obj<partition_type> partition = distributeMeshV(mesh, newMesh, renumbering, sendMeshOverlap, recvMeshOverlap);
326:       if (mesh->debug()) {
327:         std::cout << "["<<mesh->commRank()<<"]: Mesh Renumbering:" << std::endl;
328:         for(typename Mesh::renumbering_type::const_iterator r_iter = renumbering.begin(); r_iter != renumbering.end(); ++r_iter) {
329:           std::cout << "["<<mesh->commRank()<<"]:   global point " << r_iter->first << " --> " << " local point " << r_iter->second << std::endl;
330:         }
331:       }
332:       // Distribute the coordinates
333:       PETSc::Log::Event("DistributeCoords").begin();
334:       const Obj<typename Mesh::real_section_type>& coordinates         = mesh->getRealSection("coordinates");
335:       const Obj<typename Mesh::real_section_type>& parallelCoordinates = newMesh->getRealSection("coordinates");

337:       newMesh->setupCoordinates(parallelCoordinates);
338:       distributeSection(coordinates, partition, renumbering, sendMeshOverlap, recvMeshOverlap, parallelCoordinates);
339:       PETSc::Log::Event("DistributeCoords").end();
340:       // Distribute other sections
341:       if (mesh->getRealSections()->size() > 1) {
342:         PETSc::Log::Event("DistributeRealSec").begin();
343:         Obj<std::set<std::string> > names = mesh->getRealSections();
344:         int                         n     = 0;

346:         for(std::set<std::string>::const_iterator n_iter = names->begin(); n_iter != names->end(); ++n_iter) {
347:           if (*n_iter == "coordinates")   continue;
348:           std::cout << "ERROR: Did not distribute real section " << *n_iter << std::endl;
349:           ++n;
350:         }
351:         PETSc::Log::Event("DistributeRealSec").end();
352:         if (n) {throw ALE::Exception("Need to distribute more real sections");}
353:       }
354:       if (mesh->getIntSections()->size() > 0) {
355:         PETSc::Log::Event("DistributeIntSec").begin();
356:         Obj<std::set<std::string> > names = mesh->getIntSections();

358:         for(std::set<std::string>::const_iterator n_iter = names->begin(); n_iter != names->end(); ++n_iter) {
359:           const Obj<typename Mesh::int_section_type>& section    = mesh->getIntSection(*n_iter);
360:           const Obj<typename Mesh::int_section_type>& newSection = newMesh->getIntSection(*n_iter);
361: 
362:           // We assume all integer sections are complete sections
363:           newSection->setChart(newMesh->getSieve()->getChart());
364:           distributeSection(section, partition, renumbering, sendMeshOverlap, recvMeshOverlap, newSection);
365:         }
366:         PETSc::Log::Event("DistributeIntSec").end();
367:       }
368:       if (mesh->getArrowSections()->size() > 1) {
369:         throw ALE::Exception("Need to distribute more arrow sections");
370:       }
371:       // Distribute labels
372:       PETSc::Log::Event("DistributeLabels").begin();
373:       const typename Mesh::labels_type& labels = mesh->getLabels();

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

380: #ifdef IMESH_NEW_LABELS
381:         newLabel->setChart(newMesh->getSieve()->getChart());
382:         // Size the local mesh
383:         Partitioner::sizeLocalSieveV(origLabel, partition, renumbering, newLabel);
384:         // Create the remote meshes
385:         completeConesV(origLabel, newLabel, renumbering, sendMeshOverlap, recvMeshOverlap);
386:         // Create the local mesh
387:         Partitioner::createLocalSieveV(origLabel, partition, renumbering, newLabel);
388:         newLabel->symmetrize();
389: #else
390:         distributeLabelV(newMesh->getSieve(), origLabel, partition, renumbering, sendMeshOverlap, recvMeshOverlap, newLabel);
391: #endif
392:       }
393:       PETSc::Log::Event("DistributeLabels").end();
394:       // Create the parallel overlap
395:       PETSc::Log::Event("CreateOverlap").begin();
396:       Obj<typename Mesh::send_overlap_type> sendParallelMeshOverlap = newMesh->getSendOverlap();
397:       Obj<typename Mesh::recv_overlap_type> recvParallelMeshOverlap = newMesh->getRecvOverlap();
398:       //   Can I figure this out in a nicer way?
399:       ALE::SetFromMap<std::map<point_type,point_type> > globalPoints(renumbering);

401:       ALE::OverlapBuilder<>::constructOverlap(globalPoints, renumbering, sendParallelMeshOverlap, recvParallelMeshOverlap);
402:       newMesh->setCalculatedOverlap(true);
403:       PETSc::Log::Event("CreateOverlap").end();
404:     };
405:     template<typename Label, typename Partition, typename Renumbering, typename SendOverlap, typename RecvOverlap, typename NewLabel>
406:     static void distributeLabel(const Obj<typename Mesh::sieve_type>& sieve, const Obj<Label>& l, const Obj<Partition>& partition, Renumbering& renumbering, const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<NewLabel>& newL) {
407:       Partitioner::createLocalSifter(l, partition, renumbering, newL);
408:       //completeCones(l, newL, renumbering, sendMeshOverlap, recvMeshOverlap);
409:       {
410:         typedef ALE::UniformSection<point_type, int>                cones_type;
411:         typedef ALE::LabelSection<typename Mesh::sieve_type, Label> cones_wrapper_type;
412:         Obj<cones_wrapper_type> cones        = new cones_wrapper_type(sieve, l);
413:         Obj<cones_type>         overlapCones = new cones_type(l->comm(), l->debug());

415:         ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, cones, overlapCones);
416:         if (l->debug()) {overlapCones->view("Overlap Label Values");}
417:         // Inserts cones into newL (must renumber here)
418:         //ALE::Pullback::InsertionBinaryFusion::fuse(overlapCones, recvOverlap, renumbering, newSieve);
419:         {
420:           typedef typename cones_type::point_type overlap_point_type;
421:           const Obj<typename RecvOverlap::traits::baseSequence>      rPoints = recvOverlap->base();
422:           const typename RecvOverlap::traits::baseSequence::iterator rEnd    = rPoints->end();

424:           for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
425:             const Obj<typename RecvOverlap::coneSequence>& points       = recvOverlap->cone(*p_iter);
426:             const typename RecvOverlap::target_type&       localPoint   = *p_iter;
427:             const typename cones_type::point_type&         remotePoint  = points->begin().color();
428:             const overlap_point_type                       overlapPoint = overlap_point_type(remotePoint.second, remotePoint.first);
429:             const int                                      size         = overlapCones->getFiberDimension(overlapPoint);
430:             const typename cones_type::value_type         *values       = overlapCones->restrictPoint(overlapPoint);

432:             newL->clearCone(localPoint);
433:             for(int i = 0; i < size; ++i) {newL->addCone(values[i], localPoint);}
434:           }
435:         }
436:       }
437:     };
438:     template<typename Label, typename Partition, typename Renumbering, typename SendOverlap, typename RecvOverlap, typename NewLabel>
439:     static void distributeLabelV(const Obj<typename Mesh::sieve_type>& sieve, const Obj<Label>& l, const Obj<Partition>& partition, Renumbering& renumbering, const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<NewLabel>& newL) {
440:       Partitioner::createLocalSifter(l, partition, renumbering, newL);
441:       //completeCones(l, newL, renumbering, sendMeshOverlap, recvMeshOverlap);
442:       {
443:         typedef typename Label::alloc_type::template rebind<int>::other alloc_type;
444:         typedef LabelBaseSectionV<typename Mesh::sieve_type, Label, alloc_type> atlas_type;
445:         typedef ALE::UniformSection<ALE::Pair<int, point_type>, int>            cones_type;
446:         typedef ALE::LabelSection<typename Mesh::sieve_type, Label, alloc_type, atlas_type> cones_wrapper_type;
447:         Obj<cones_wrapper_type> cones        = new cones_wrapper_type(sieve, l);
448:         Obj<cones_type>         overlapCones = new cones_type(l->comm(), l->debug());

450:         ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, cones, overlapCones);
451:         if (l->debug()) {overlapCones->view("Overlap Label Values");}
452:         // Inserts cones into newL (must renumber here)
453:         //ALE::Pullback::InsertionBinaryFusion::fuse(overlapCones, recvOverlap, renumbering, newSieve);
454:         {
455:           typedef typename cones_type::point_type overlap_point_type;
456:           const Obj<typename RecvOverlap::traits::baseSequence> rPoints = recvOverlap->base();

458:           for(typename RecvOverlap::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rPoints->end(); ++p_iter) {
459:             const Obj<typename RecvOverlap::coneSequence>& points       = recvOverlap->cone(*p_iter);
460:             const typename RecvOverlap::target_type&       localPoint   = *p_iter;
461:             const typename cones_type::point_type&         remotePoint  = points->begin().color();
462:             const overlap_point_type                       overlapPoint = overlap_point_type(remotePoint.second, remotePoint.first);
463:             const int                                      size         = overlapCones->getFiberDimension(overlapPoint);
464:             const typename cones_type::value_type         *values       = overlapCones->restrictPoint(overlapPoint);

466:             newL->clearCone(localPoint);
467:             for(int i = 0; i < size; ++i) {newL->addCone(values[i], localPoint);}
468:           }
469:         }
470:       }
471:     };
472:     template<typename Section, typename Partition, typename Renumbering, typename SendOverlap, typename RecvOverlap, typename NewSection>
473:     static void distributeSection(const Obj<Section>& s, const Obj<Partition>& partition, Renumbering& renumbering, const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<NewSection>& newS) {
474:       Partitioner::createLocalSection(s, partition, renumbering, newS);
475:       ALE::Completion::completeSection(sendOverlap, recvOverlap, s, newS);
476:     };
477:     template<typename NewMesh, typename Renumbering, typename SendOverlap, typename RecvOverlap>
478:     static Obj<partition_type> unifyMesh(const Obj<Mesh>& mesh, const Obj<NewMesh>& newMesh, Renumbering& renumbering, const Obj<SendOverlap>& sendMeshOverlap, const Obj<RecvOverlap>& recvMeshOverlap) {
479:       const Obj<partition_type> cellPartition = new partition_type(mesh->comm(), 0, mesh->commSize(), mesh->debug());
480:       const Obj<partition_type> partition     = new partition_type(mesh->comm(), 0, mesh->commSize(), mesh->debug());
481:       const Obj<typename Mesh::label_sequence>&     cells  = mesh->heightStratum(0);
482:       const typename Mesh::label_sequence::iterator cEnd   = cells->end();
483:       typename Mesh::point_type                    *values = new typename Mesh::point_type[cells->size()];
484:       int                                           c      = 0;

486:       cellPartition->setFiberDimension(0, cells->size());
487:       cellPartition->allocatePoint();
488:       for(typename Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != cEnd; ++c_iter, ++c) {
489:         values[c] = *c_iter;
490:       }
491:       cellPartition->updatePoint(0, values);
492:       delete [] values;
493:       // Close the partition over sieve points
494:       Partitioner::createPartitionClosure(mesh, cellPartition, partition);
495:       // Create the remote meshes
496:       completeMesh(mesh, partition, renumbering, newMesh, sendMeshOverlap, recvMeshOverlap);
497:       // Create the local mesh
498:       Partitioner::createLocalMesh(mesh, partition, renumbering, newMesh);
499:       newMesh->stratify();
500:       newMesh->view("Unified mesh");
501:       return partition;
502:     };
503:     static Obj<Mesh> unifyMesh(const Obj<Mesh>& mesh) {
504:       typedef ALE::Sifter<point_type,rank_type,point_type> mesh_send_overlap_type;
505:       typedef ALE::Sifter<rank_type,point_type,point_type> mesh_recv_overlap_type;
506:       const Obj<Mesh>                      newMesh         = new Mesh(mesh->comm(), mesh->getDimension(), mesh->debug());
507:       const Obj<typename Mesh::sieve_type> newSieve        = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
508:       const Obj<mesh_send_overlap_type>    sendMeshOverlap = new mesh_send_overlap_type(mesh->comm(), mesh->debug());
509:       const Obj<mesh_recv_overlap_type>    recvMeshOverlap = new mesh_recv_overlap_type(mesh->comm(), mesh->debug());
510:       std::map<point_type,point_type>      renumbering;

512:       newMesh->setSieve(newSieve);
513:       const Obj<partition_type> partition = unifyMesh(mesh, newMesh, renumbering, sendMeshOverlap, recvMeshOverlap);
514:       // Unify coordinates
515:       const Obj<typename Mesh::real_section_type>& coordinates    = mesh->getRealSection("coordinates");
516:       const Obj<typename Mesh::real_section_type>& newCoordinates = newMesh->getRealSection("coordinates");

518:       newMesh->setupCoordinates(newCoordinates);
519:       distributeSection(coordinates, partition, renumbering, sendMeshOverlap, recvMeshOverlap, newCoordinates);
520:       // Unify labels
521:       const typename Mesh::labels_type& labels = mesh->getLabels();

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

528:         //completeCones(label, newLabel, renumbering, sendMeshOverlap, recvMeshOverlap);
529:         {
530:           typedef ALE::UniformSection<point_type, int> cones_type;
531:           typedef ALE::LabelSection<typename Mesh::sieve_type,typename Mesh::label_type> cones_wrapper_type;
532:           Obj<cones_wrapper_type> cones        = new cones_wrapper_type(mesh->getSieve(), label);
533:           Obj<cones_type>         overlapCones = new cones_type(label->comm(), label->debug());

535:           ALE::Pullback::SimpleCopy::copy(sendMeshOverlap, recvMeshOverlap, cones, overlapCones);
536:           if (label->debug()) {overlapCones->view("Overlap Label Values");}
537:           // Inserts cones into parallelMesh (must renumber here)
538:           //ALE::Pullback::InsertionBinaryFusion::fuse(overlapCones, recvMeshOverlap, renumbering, newSieve);
539:           {
540:             const Obj<typename mesh_recv_overlap_type::traits::baseSequence> rPoints = recvMeshOverlap->base();

542:             for(typename mesh_recv_overlap_type::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rPoints->end(); ++p_iter) {
543:               const Obj<typename mesh_recv_overlap_type::coneSequence>& points      = recvMeshOverlap->cone(*p_iter);
544:               const typename mesh_recv_overlap_type::target_type&       localPoint  = *p_iter;
545:               const typename cones_type::point_type&                    remotePoint = points->begin().color();
546:               const int                                                 size        = overlapCones->getFiberDimension(remotePoint);
547:               const typename cones_type::value_type                    *values      = overlapCones->restrictPoint(remotePoint);

549:               newLabel->clearCone(localPoint);
550:               for(int i = 0; i < size; ++i) {newLabel->addCone(values[i], localPoint);}
551:             }
552:           }
553:         }
554:         //newLabel->add(label, newSieve);
555:         Partitioner::createLocalSifter(label, partition, renumbering, newLabel);
556:         newLabel->view(l_iter->first.c_str());
557:       }
558:       return newMesh;
559:     };
560:   };
561:   template<typename Bundle_>
562:   class Distribution {
563:   public:
564:     typedef Bundle_                                                                     bundle_type;
565:     typedef typename bundle_type::sieve_type                                            sieve_type;
566:     typedef typename bundle_type::point_type                                            point_type;
567:     typedef typename bundle_type::alloc_type                                            alloc_type;
568:     typedef typename bundle_type::send_overlap_type                                     send_overlap_type;
569:     typedef typename bundle_type::recv_overlap_type                                     recv_overlap_type;
570:     typedef typename ALE::New::Completion<bundle_type, typename sieve_type::point_type>                            sieveCompletion;
571:     typedef typename ALE::New::SectionCompletion<bundle_type, typename bundle_type::real_section_type::value_type> sectionCompletion;
572:   public:
575:     static void createPartitionOverlap(const Obj<bundle_type>& bundle, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
576:       const Obj<send_overlap_type>& topSendOverlap = bundle->getSendOverlap();
577:       const Obj<recv_overlap_type>& topRecvOverlap = bundle->getRecvOverlap();
578:       const Obj<typename send_overlap_type::traits::baseSequence> base = topSendOverlap->base();
579:       const Obj<typename recv_overlap_type::traits::capSequence>  cap  = topRecvOverlap->cap();
580:       const int rank = bundle->commRank();

582:       if (base->empty()) {
583:         if (rank == 0) {
584:           for(int p = 1; p < bundle->commSize(); p++) {
585:             // The arrow is from local partition point p (source) to remote partition point p (color) on rank p (target)
586:             sendOverlap->addCone(p, p, p);
587:           }
588:         }
589:       } else {
590:         for(typename send_overlap_type::traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
591:           const int& p = *b_iter;
592:           // The arrow is from local partition point p (source) to remote partition point p (color) on rank p (target)
593:           sendOverlap->addCone(p, p, p);
594:         }
595:       }
596:       if (cap->empty()) {
597:         if (rank != 0) {
598:           // The arrow is from local partition point rank (color) on rank 0 (source) to remote partition point rank (target)
599:           recvOverlap->addCone(0, rank, rank);
600:         }
601:       } else {
602:         for(typename recv_overlap_type::traits::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
603:           const int& p = *c_iter;
604:           // The arrow is from local partition point rank (color) on rank p (source) to remote partition point rank (target)
605:           recvOverlap->addCone(p, rank, rank);
606:         }
607:       }
608:     };
611:     template<typename Partitioner>
612:     static typename Partitioner::part_type *createAssignment(const Obj<bundle_type>& bundle, const int dim, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const int height = 0) {
613:       // 1) Form partition point overlap a priori
614:       createPartitionOverlap(bundle, sendOverlap, recvOverlap);
615:       if (bundle->debug()) {
616:         sendOverlap->view("Send overlap for partition");
617:         recvOverlap->view("Receive overlap for partition");
618:       }
619:       // 2) Partition the mesh
620:       if (height == 0) {
621:         return Partitioner::partitionSieve(bundle, dim);
622:       } else if (height == 1) {
623:         return Partitioner::partitionSieveByFace(bundle, dim);
624:       }
625:       throw ALE::Exception("Invalid partition height");
626:     };
629:     // Partition a bundle on process 0 and scatter to all processes
630:     static void scatterBundle(const Obj<bundle_type>& bundle, const int dim, const Obj<bundle_type>& bundleNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const std::string& partitioner, const int height = 0, const Obj<bundle_type>& subBundle = NULL, const Obj<bundle_type>& subBundleNew = NULL) {
631:       if (height == 0) {
632:         if (partitioner == "chaco") {
633: #ifdef PETSC_HAVE_CHACO
634:           typedef typename ALE::New::Chaco::Partitioner<bundle_type> Partitioner;
635:           typedef typename ALE::New::Partitioner<bundle_type>        GenPartitioner;
636:           typedef typename Partitioner::part_type                    part_type;

638:           part_type *assignment = scatterBundle<Partitioner>(bundle, dim, bundleNew, sendOverlap, recvOverlap, height);
639:           if (!subBundle.isNull() && !subBundleNew.isNull()) {
640:             part_type *subAssignment = GenPartitioner::subordinatePartition(bundle, 1, subBundle, assignment);
641:             const Obj<sieve_type>& sieve      = subBundle->getSieve();
642:             const Obj<sieve_type>& sieveNew   = new Mesh::sieve_type(subBundle->comm(), subBundle->debug());
643:             const int              numCells   = subBundle->heightStratum(height)->size();

645:             subBundleNew->setSieve(sieveNew);
646:             sieveCompletion::scatterSieve(subBundle, sieve, dim, sieveNew, sendOverlap, recvOverlap, height, numCells, subAssignment);
647:             subBundleNew->stratify();
648:             if (subAssignment != NULL) delete [] subAssignment;
649:           }
650:           if (assignment != NULL) delete [] assignment;
651: #else
652:           throw ALE::Exception("Chaco is not installed. Reconfigure with the flag --download-chaco");
653: #endif
654:         } else if (partitioner == "parmetis") {
655: #ifdef PETSC_HAVE_PARMETIS
656:           typedef typename ALE::New::ParMetis::Partitioner<bundle_type> Partitioner;
657:           typedef typename ALE::New::Partitioner<bundle_type>           GenPartitioner;
658:           typedef typename Partitioner::part_type                       part_type;

660:           part_type *assignment = scatterBundle<Partitioner>(bundle, dim, bundleNew, sendOverlap, recvOverlap, height);
661:           if (!subBundle.isNull() && !subBundleNew.isNull()) {
662:             part_type *subAssignment = GenPartitioner::subordinatePartition(bundle, 1, subBundle, assignment);
663:             const Obj<sieve_type>& sieve      = subBundle->getSieve();
664:             const Obj<sieve_type>& sieveNew   = new Mesh::sieve_type(subBundle->comm(), subBundle->debug());
665:             const int              numCells   = subBundle->heightStratum(height)->size();

667:             subBundleNew->setSieve(sieveNew);
668:             sieveCompletion::scatterSieve(subBundle, sieve, dim, sieveNew, sendOverlap, recvOverlap, height, numCells, subAssignment);
669:             subBundleNew->stratify();
670:             if (subAssignment != NULL) delete [] subAssignment;
671:           }
672:           if (assignment != NULL) delete [] assignment;
673: #else
674:           throw ALE::Exception("ParMetis is not installed. Reconfigure with the flag --download-parmetis");
675: #endif
676:         } else {
677:           throw ALE::Exception("Unknown partitioner");
678:         }
679:       } else if (height == 1) {
680:         if (partitioner == "zoltan") {
681: #ifdef PETSC_HAVE_ZOLTAN
682:           typedef typename ALE::New::Zoltan::Partitioner<bundle_type> Partitioner;
683:           typedef typename Partitioner::part_type                     part_type;

685:           part_type *assignment = scatterBundle<Partitioner>(bundle, dim, bundleNew, sendOverlap, recvOverlap, height);
686:           if (assignment != NULL) delete [] assignment;
687: #else
688:           throw ALE::Exception("Zoltan is not installed. Reconfigure with the flag --download-zoltan");
689: #endif
690:         } else if (partitioner == "parmetis") {
691: #ifdef PETSC_HAVE_PARMETIS
692:           typedef typename ALE::New::ParMetis::Partitioner<bundle_type> Partitioner;
693:           typedef typename Partitioner::part_type                       part_type;

695:           part_type *assignment = scatterBundle<Partitioner>(bundle, dim, bundleNew, sendOverlap, recvOverlap, height);
696:           if (assignment != NULL) delete [] assignment;
697: #else
698:           throw ALE::Exception("ParMetis is not installed. Reconfigure with the flag --download-parmetis");
699: #endif
700:         } else {
701:           throw ALE::Exception("Unknown partitioner");
702:         }
703:       } else {
704:         throw ALE::Exception("Invalid partition height");
705:       }
706:     };
707:     template<typename Partitioner>
708:     static typename Partitioner::part_type *scatterBundle(const Obj<bundle_type>& bundle, const int dim, const Obj<bundle_type>& bundleNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const int height = 0) {
709:       typename Partitioner::part_type *assignment = createAssignment<Partitioner>(bundle, dim, sendOverlap, recvOverlap, height);
710:       const Obj<sieve_type>&           sieve      = bundle->getSieve();
711:       const Obj<sieve_type>&           sieveNew   = bundleNew->getSieve();
712:       const int                        numPoints  = bundle->heightStratum(height)->size();

714:       sieveCompletion::scatterSieve(bundle, sieve, dim, sieveNew, sendOverlap, recvOverlap, height, numPoints, assignment);
715:       bundleNew->stratify();
716:       return assignment;
717:     };
720:     static Obj<ALE::Mesh> distributeMesh(const Obj<ALE::Mesh>& serialMesh, const int height = 0, const std::string& partitioner = "chaco") {
721:       MPI_Comm                          comm          = serialMesh->comm();
722:       const int                         dim           = serialMesh->getDimension();
723:       Obj<ALE::Mesh>                    parallelMesh  = new ALE::Mesh(comm, dim, serialMesh->debug());
724:       const Obj<ALE::Mesh::sieve_type>& parallelSieve = new ALE::Mesh::sieve_type(comm, serialMesh->debug());

726:       ALE_LOG_EVENT_BEGIN;
727:       parallelMesh->setSieve(parallelSieve);
728:       if (serialMesh->debug()) {serialMesh->view("Serial mesh");}

730:       // Distribute cones
731:       Obj<send_overlap_type> sendOverlap = new send_overlap_type(comm, serialMesh->debug());
732:       Obj<recv_overlap_type> recvOverlap = new recv_overlap_type(comm, serialMesh->debug());
733:       scatterBundle(serialMesh, dim, parallelMesh, sendOverlap, recvOverlap, partitioner, height);
734:       parallelMesh->setDistSendOverlap(sendOverlap);
735:       parallelMesh->setDistRecvOverlap(recvOverlap);

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

740:       for(typename bundle_type::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
741:         if (parallelMesh->hasLabel(l_iter->first)) continue;
742:         const Obj<typename bundle_type::label_type>& serialLabel   = l_iter->second;
743:         const Obj<typename bundle_type::label_type>& parallelLabel = parallelMesh->createLabel(l_iter->first);
744:         // Create local label
745: #define NEW_LABEL
746: #ifdef NEW_LABEL
747:         parallelLabel->add(serialLabel, parallelSieve);
748: #else
749:         const Obj<typename bundle_type::label_type::traits::baseSequence>& base = serialLabel->base();

751:         for(typename bundle_type::label_type::traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
752:           if (parallelSieve->capContains(*b_iter) || parallelSieve->baseContains(*b_iter)) {
753:             parallelLabel->addArrow(*serialLabel->cone(*b_iter)->begin(), *b_iter);
754:           }
755:         }
756: #endif
757:         // Get remote labels
758:         sieveCompletion::scatterCones(serialLabel, parallelLabel, sendOverlap, recvOverlap);
759:       }

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

764:       for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
765:         parallelMesh->setRealSection(*name, distributeSection(serialMesh->getRealSection(*name), parallelMesh, sendOverlap, recvOverlap));
766:       }
767:       sections = serialMesh->getIntSections();
768:       for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
769:         parallelMesh->setIntSection(*name, distributeSection(serialMesh->getIntSection(*name), parallelMesh, sendOverlap, recvOverlap));
770:       }
771:       sections = serialMesh->getArrowSections();

773:       for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
774:         parallelMesh->setArrowSection(*name, distributeArrowSection(serialMesh->getArrowSection(*name), serialMesh, parallelMesh, sendOverlap, recvOverlap));
775:       }
776:       if (parallelMesh->debug()) {parallelMesh->view("Parallel Mesh");}
777:       ALE_LOG_EVENT_END;
778:       return parallelMesh;
779:     };
782:     template<typename Section>
783:     static void updateSectionLocal(const Obj<Section>& oldSection, const Obj<bundle_type>& newBundle, const Obj<Section>& newSection) {
784:       const Obj<typename bundle_type::sieve_type>&    newSieve = newBundle->getSieve();
785:       const typename Section::atlas_type::chart_type& oldChart = oldSection->getChart();

787:       for(typename Section::atlas_type::chart_type::const_iterator c_iter = oldChart.begin(); c_iter != oldChart.end(); ++c_iter) {
788:         if (newSieve->capContains(*c_iter) || newSieve->baseContains(*c_iter)) {
789:           newSection->setFiberDimension(*c_iter, oldSection->getFiberDimension(*c_iter));
790:         }
791:       }
792:       newBundle->allocate(newSection);
793:       const typename Section::atlas_type::chart_type& newChart = newSection->getChart();

795:       for(typename Section::atlas_type::chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
796:         newSection->updatePointAll(*c_iter, oldSection->restrictPoint(*c_iter));
797:       }
798:     };
801:     template<typename RecvSection, typename Section>
802:     static void updateSectionRemote(const Obj<recv_overlap_type>& recvOverlap, const Obj<RecvSection>& recvSection, const Obj<bundle_type>& newBundle, const Obj<Section>& newSection) {
803:       Obj<typename recv_overlap_type::traits::baseSequence> recvPoints = recvOverlap->base();

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

809:         for(typename recv_overlap_type::traits::coneSequence::iterator p_iter = recvPatches->begin(); p_iter != end; ++p_iter) {
810:           newSection->addPoint(*r_iter, recvSection->getSection(*p_iter)->getFiberDimension(*r_iter));
811:         }
812:       }
813:       newBundle->reallocate(newSection);
814:       for(typename recv_overlap_type::traits::baseSequence::iterator r_iter = recvPoints->begin(); r_iter != recvPoints->end(); ++r_iter) {
815:         const Obj<typename recv_overlap_type::traits::coneSequence>&     recvPatches = recvOverlap->cone(*r_iter);
816:         const typename recv_overlap_type::traits::coneSequence::iterator end         = recvPatches->end();

818:         for(typename recv_overlap_type::traits::coneSequence::iterator p_iter = recvPatches->begin(); p_iter != end; ++p_iter) {
819:           if (recvSection->getSection(*p_iter)->getFiberDimension(*r_iter)) {
820:             newSection->updatePointAll(*r_iter, recvSection->getSection(*p_iter)->restrictPoint(*r_iter));
821:           }
822:         }
823:       }
824:     };
827:     template<typename Section>
828:     static Obj<Section> distributeSection(const Obj<Section>& serialSection, const Obj<bundle_type>& parallelBundle, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
829:       if (serialSection->debug()) {
830:         serialSection->view("Serial Section");
831:       }
832:       typedef typename alloc_type::template rebind<typename Section::value_type>::other value_alloc_type;
833:       typedef ALE::Field<send_overlap_type, int, ALE::Section<point_type, typename Section::value_type, value_alloc_type> > send_section_type;
834:       typedef ALE::Field<recv_overlap_type, int, ALE::Section<point_type, typename Section::value_type, value_alloc_type> > recv_section_type;
835:       typedef ALE::New::SizeSection<Section> SectionSizer;
836:       Obj<Section>                 parallelSection = new Section(serialSection->comm(), serialSection->debug());
837:       const Obj<send_section_type> sendSection     = new send_section_type(serialSection->comm(), serialSection->debug());
838:       const Obj<recv_section_type> recvSection     = new recv_section_type(serialSection->comm(), sendSection->getTag(), serialSection->debug());
839:       const Obj<SectionSizer>      sizer           = new SectionSizer(serialSection);

841:       updateSectionLocal(serialSection, parallelBundle, parallelSection);
842:       sectionCompletion::completeSection(sendOverlap, recvOverlap, sizer, serialSection, sendSection, recvSection);
843:       updateSectionRemote(recvOverlap, recvSection, parallelBundle, parallelSection);
844:       if (parallelSection->debug()) {
845:         parallelSection->view("Parallel Section");
846:       }
847:       return parallelSection;
848:     };
851:     template<typename Section>
852:     static void updateArrowSectionLocal(const Obj<Section>& oldSection, const Obj<bundle_type>& newBundle, const Obj<Section>& newSection) {
853:       const Obj<typename bundle_type::sieve_type>&    newSieve = newBundle->getSieve();
854:       const typename Section::atlas_type::chart_type& oldChart = oldSection->getChart();

856:       for(typename Section::atlas_type::chart_type::const_iterator c_iter = oldChart.begin(); c_iter != oldChart.end(); ++c_iter) {
857:         // Dmitry should provide a Sieve::contains(MinimalArrow) method
858:         if (newSieve->capContains(c_iter->source) && newSieve->baseContains(c_iter->target)) {
859:           newSection->setFiberDimension(*c_iter, oldSection->getFiberDimension(*c_iter));
860:         }
861:       }
862:       //newBundle->allocate(newSection);
863:       const typename Section::atlas_type::chart_type& newChart = newSection->getChart();

865:       for(typename Section::atlas_type::chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
866:         newSection->updatePointAll(*c_iter, oldSection->restrictPoint(*c_iter));
867:       }
868:     };
871:     template<typename RecvSection, typename Section>
872:     static void updateArrowSectionRemote(const Obj<recv_overlap_type>& recvOverlap, const Obj<RecvSection>& recvSection, const Obj<bundle_type>& newBundle, const Obj<Section>& newSection) {
873:       Obj<typename recv_overlap_type::traits::baseSequence> recvPoints = recvOverlap->base();

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

879:         for(typename bundle_type::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
880:           newSection->setFiberDimension(typename Section::point_type(*c_iter, *r_iter), 1);
881:         }
882:       }
883:       //newBundle->reallocate(newSection);
884:       for(typename recv_overlap_type::traits::baseSequence::iterator r_iter = recvPoints->begin(); r_iter != recvPoints->end(); ++r_iter) {
885:         const Obj<typename recv_overlap_type::traits::coneSequence>&     recvPatches = recvOverlap->cone(*r_iter);
886:         const typename recv_overlap_type::traits::coneSequence::iterator recvEnd     = recvPatches->end();

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

891:           if (section->getFiberDimension(*r_iter)) {
892:             const Obj<typename bundle_type::sieve_type::traits::coneSequence>&     cone    = newBundle->getSieve()->cone(*r_iter);
893:             const typename bundle_type::sieve_type::traits::coneSequence::iterator end     = cone->end();
894:             const typename RecvSection::value_type                                *values  = section->restrictPoint(*r_iter);
895:             int                                                                    c       = -1;

897:             for(typename bundle_type::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
898:               newSection->updatePoint(typename Section::point_type(*c_iter, *r_iter), &values[++c]);
899:             }
900:           }
901:         }
902:       }
903:     };
906:     template<typename Section>
907:     static Obj<Section> distributeArrowSection(const Obj<Section>& serialSection, const Obj<bundle_type>& serialBundle, const Obj<bundle_type>& parallelBundle, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
908:       if (serialSection->debug()) {
909:         serialSection->view("Serial ArrowSection");
910:       }
911:       typedef typename alloc_type::template rebind<typename Section::value_type>::other value_alloc_type;
912:       typedef ALE::Field<send_overlap_type, int, ALE::Section<point_type, typename Section::value_type, value_alloc_type> > send_section_type;
913:       typedef ALE::Field<recv_overlap_type, int, ALE::Section<point_type, typename Section::value_type, value_alloc_type> > recv_section_type;
914:       typedef ALE::New::ConeSizeSection<bundle_type, sieve_type> SectionSizer;
915:       typedef ALE::New::ArrowSection<sieve_type, Section>        ArrowFiller;
916:       Obj<Section>                 parallelSection = new Section(serialSection->comm(), serialSection->debug());
917:       const Obj<send_section_type> sendSection     = new send_section_type(serialSection->comm(), serialSection->debug());
918:       const Obj<recv_section_type> recvSection     = new recv_section_type(serialSection->comm(), sendSection->getTag(), serialSection->debug());
919:       const Obj<SectionSizer>      sizer           = new SectionSizer(serialBundle, serialBundle->getSieve());
920:       const Obj<ArrowFiller>       filler          = new ArrowFiller(serialBundle->getSieve(), serialSection);

922:       updateArrowSectionLocal(serialSection, parallelBundle, parallelSection);
923:       sectionCompletion::completeSection(sendOverlap, recvOverlap, sizer, filler, sendSection, recvSection);
924:       updateArrowSectionRemote(recvOverlap, recvSection, parallelBundle, parallelSection);
925:       if (parallelSection->debug()) {
926:         parallelSection->view("Parallel ArrowSection");
927:       }
928:       return parallelSection;
929:     };
930:     static void unifyBundle(const Obj<bundle_type>& bundle, const int dim, const Obj<bundle_type>& bundleNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
931:       typedef int part_type;
932:       const Obj<sieve_type>& sieve    = bundle->getSieve();
933:       const Obj<sieve_type>& sieveNew = bundleNew->getSieve();
934:       const int              rank     = bundle->commRank();
935:       const int              debug    = bundle->debug();

937:       // 1) Form partition point overlap a priori
938:       if (rank == 0) {
939:         for(int p = 1; p < sieve->commSize(); p++) {
940:           // The arrow is from remote partition point 0 on rank p to local partition point 0
941:           recvOverlap->addCone(p, 0, 0);
942:         }
943:       } else {
944:         // The arrow is from local partition point 0 to remote partition point 0 on rank 0
945:         sendOverlap->addCone(0, 0, 0);
946:       }
947:       if (debug) {
948:         sendOverlap->view("Send overlap for partition");
949:         recvOverlap->view("Receive overlap for partition");
950:       }
951:       // 2) Partition the mesh
952:       int        numCells   = bundle->heightStratum(0)->size();
953:       part_type *assignment = new part_type[numCells];

955:       for(int c = 0; c < numCells; ++c) {
956:         assignment[c] = 0;
957:       }
958:       // 3) Scatter the sieve
959:       sieveCompletion::scatterSieve(bundle, sieve, dim, sieveNew, sendOverlap, recvOverlap, 0, numCells, assignment);
960:       bundleNew->stratify();
961:       // 4) Cleanup
962:       if (assignment != NULL) delete [] assignment;
963:     };
966:     static Obj<ALE::Mesh> unifyMesh(const Obj<ALE::Mesh>& parallelMesh) {
967:       const int                         dim         = parallelMesh->getDimension();
968:       Obj<ALE::Mesh>                    serialMesh  = new ALE::Mesh(parallelMesh->comm(), dim, parallelMesh->debug());
969:       const Obj<ALE::Mesh::sieve_type>& serialSieve = new ALE::Mesh::sieve_type(parallelMesh->comm(), parallelMesh->debug());

971:       ALE_LOG_EVENT_BEGIN;
972:       serialMesh->setSieve(serialSieve);
973:       if (parallelMesh->debug()) {
974:         parallelMesh->view("Parallel topology");
975:       }

977:       // Unify cones
978:       Obj<send_overlap_type> sendOverlap = new send_overlap_type(serialMesh->comm(), serialMesh->debug());
979:       Obj<recv_overlap_type> recvOverlap = new recv_overlap_type(serialMesh->comm(), serialMesh->debug());
980:       unifyBundle(parallelMesh, dim, serialMesh, sendOverlap, recvOverlap);
981:       serialMesh->setDistSendOverlap(sendOverlap);
982:       serialMesh->setDistRecvOverlap(recvOverlap);

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

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

992:         sieveCompletion::scatterCones(parallelLabel, serialLabel, sendOverlap, recvOverlap);
993:       }

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

998:       for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
999:         serialMesh->setRealSection(*name, distributeSection(parallelMesh->getRealSection(*name), serialMesh, sendOverlap, recvOverlap));
1000:       }
1001:       sections = parallelMesh->getIntSections();
1002:       for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
1003:         serialMesh->setIntSection(*name, distributeSection(parallelMesh->getIntSection(*name), serialMesh, sendOverlap, recvOverlap));
1004:       }
1005:       sections = parallelMesh->getArrowSections();
1006:       for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
1007:         serialMesh->setArrowSection(*name, distributeArrowSection(parallelMesh->getArrowSection(*name), parallelMesh, serialMesh, sendOverlap, recvOverlap));
1008:       }
1009:       if (serialMesh->debug()) {serialMesh->view("Serial Mesh");}
1010:       ALE_LOG_EVENT_END;
1011:       return serialMesh;
1012:     };
1013:   public: // Do not like these
1016:     // This is just crappy. We could introduce another phase to find out exactly what
1017:     //   indices people do not have in the global order after communication
1018:     template<typename OrigSendOverlap, typename OrigRecvOverlap, typename SendSection, typename RecvSection>
1019:     static void updateOverlap(const Obj<OrigSendOverlap>& origSendOverlap, const Obj<OrigRecvOverlap>& origRecvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
1020:       const typename SendSection::sheaf_type& sendRanks = sendSection->getPatches();
1021:       const typename RecvSection::sheaf_type& recvRanks = recvSection->getPatches();

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

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

1032:           for(int p = 0; p < size; p++) {
1033:             if (origSendOverlap->support(points[p])->size() == 0) {
1034:               sendOverlap->addArrow(points[p], rank, points[p]);
1035:             }
1036:           }
1037:         }
1038:       }
1039:       for(typename RecvSection::sheaf_type::const_iterator p_iter = recvRanks.begin(); p_iter != recvRanks.end(); ++p_iter) {
1040:         const typename RecvSection::patch_type&               rank    = p_iter->first;
1041:         const Obj<typename RecvSection::section_type>&        section = p_iter->second;
1042:         const typename RecvSection::section_type::chart_type& chart   = section->getChart();

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

1048:           for(int p = 0; p < size; p++) {
1049:             if (origRecvOverlap->support(rank, points[p])->size() == 0) {
1050:               recvOverlap->addArrow(rank, points[p], points[p]);
1051:             }
1052:           }
1053:         }
1054:       }
1055:     };
1058:     template<typename RecvOverlap, typename RecvSection>
1059:     static void updateSieve(const Obj<RecvOverlap>& recvOverlap, const Obj<RecvSection>& recvSection, const Obj<sieve_type>& sieve) {
1060: #if 1
1061:       Obj<typename RecvOverlap::traits::baseSequence> recvPoints = recvOverlap->base();

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

1067:         for(typename RecvOverlap::traits::coneSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
1068:           const typename RecvOverlap::target_type&       remotePoint = r_iter.color();
1069:           const int                                      rank        = *r_iter;
1070:           const Obj<typename RecvSection::section_type>& section     = recvSection->getSection(rank);
1071:           const typename RecvSection::value_type        *points      = section->restrictPoint(remotePoint);
1072:           const int                                      size        = section->getFiberDimension(remotePoint);
1073:           int                                            c           = 0;

1075:           ///std::cout << "["<<recvSection->commRank()<<"]: Receiving " << size << " points from rank " << rank << std::endl;
1076:           for(int p = 0; p < size; p++) {
1077:             // rank -- remote point --> local point
1078:             if (recvOverlap->support(rank, points[p])->size()) {
1079:               sieve->addArrow(*recvOverlap->support(rank, points[p])->begin(), localPoint, c);
1080:               ///std::cout << "["<<recvSection->commRank()<<"]:   1Adding arrow " << *recvOverlap->support(rank, points[p])->begin() << "("<<points[p]<<") --> " << localPoint << std::endl;
1081:             } else {
1082:               sieve->addArrow(points[p], localPoint, c);
1083:               ///std::cout << "["<<recvSection->commRank()<<"]:   2Adding arrow " << points[p] << " --> " << localPoint << std::endl;
1084:             }
1085:           }
1086:         }
1087:       }
1088: #else
1089:       const typename RecvSection::sheaf_type& ranks = recvSection->getPatches();

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

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

1100:           std::cout << "["<<recvSection->commRank()<<"]: Receiving " << size << " points from rank " << p_iter->first << std::endl;
1101:           for(int p = 0; p < size; p++) {
1102:             //sieve->addArrow(points[p], *b_iter, c++);
1103:             sieve->addArrow(points[p], *b_iter, c);
1104:             std::cout << "["<<recvSection->commRank()<<"]:   Adding arrow " << points[p] << " --> " << *b_iter << std::endl;
1105:           }
1106:         }
1107:       }
1108: #endif
1109:     };
1112:     template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
1113:     static void coneCompletion(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<bundle_type>& bundle, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
1114:       if (sendOverlap->commSize() == 1) return;
1115:       // Distribute cones
1116:       const Obj<sieve_type>&                                 sieve           = bundle->getSieve();
1117:       const Obj<typename sieveCompletion::topology_type>     secTopology     = sieveCompletion::completion::createSendTopology(sendOverlap);
1118:       const Obj<typename sieveCompletion::cone_size_section> coneSizeSection = new typename sieveCompletion::cone_size_section(bundle, sieve);
1119:       const Obj<typename sieveCompletion::cone_section>      coneSection     = new typename sieveCompletion::cone_section(sieve);
1120:       sieveCompletion::completion::completeSection(sendOverlap, recvOverlap, coneSizeSection, coneSection, sendSection, recvSection);
1121:       // Update cones
1122:       updateSieve(recvOverlap, recvSection, sieve);
1123:     };
1126:     template<typename Section>
1127:     static void completeSection(const Obj<bundle_type>& bundle, const Obj<Section>& section) {
1128:       typedef typename Distribution<bundle_type>::sieveCompletion sieveCompletion;
1129:       typedef typename bundle_type::send_overlap_type             send_overlap_type;
1130:       typedef typename bundle_type::recv_overlap_type             recv_overlap_type;
1131:       typedef typename Section::value_type                        value_type;
1132:       typedef typename alloc_type::template rebind<typename Section::value_type>::other value_alloc_type;
1133:       typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
1134:       typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
1135:       typedef ALE::New::SizeSection<Section>                                SectionSizer;
1136:       const int debug = section->debug();

1138:       bundle->constructOverlap();
1139:       const Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
1140:       const Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();
1141:       const Obj<send_section_type> sendSection = new send_section_type(section->comm(), section->debug());
1142:       const Obj<recv_section_type> recvSection = new recv_section_type(section->comm(), sendSection->getTag(), section->debug());
1143:       const Obj<SectionSizer>      sizer       = new SectionSizer(section);

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

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

1153:         for(typename recv_overlap_type::traits::coneSequence::iterator p_iter = recvPatches->begin(); p_iter != end; ++p_iter) {
1154:           if (recvSection->getSection(*p_iter)->getFiberDimension(p_iter.color())) {
1155:             if (debug) {std::cout << "["<<section->commRank()<<"]Completed point " << *r_iter << std::endl;}
1156:             section->updateAddPoint(*r_iter, recvSection->getSection(*p_iter)->restrictPoint(p_iter.color()));
1157:           }
1158:         }
1159:       }
1160:     };
1161:   };
1162: }
1163: #endif