Actual source code: Distribution.hh
petsc-3.4.2 2013-07-02

      1: #ifndef included_ALE_Distribution_hh
      2: #define included_ALE_Distribution_hh

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

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

     12: #ifndef  included_ALE_SectionCompletion_hh
     13: #include <sieve/SectionCompletion.hh>
     14: #endif

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

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

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

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

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

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

    135:       for(point_type p = 0; p < numLocalPoints; ++p) {
    136:         renumbering[localPoints[p]] = p;
    137:       }
    138:       const Obj<typename part_recv_overlap_type::traits::baseSequence> rPoints    = recvOverlap->base();
    139:       point_type                                                       localPoint = numLocalPoints;

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

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

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

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

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

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

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

    205:       // Create the cell partition
    206:       Partitioner::createPartition(mesh, cellPartition, height);
    207:       if (mesh->debug()) {
    208:         PetscViewer    viewer;

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

    232:       // Distribute the coordinates
    233:       const Obj<typename Mesh::real_section_type>& coordinates         = mesh->getRealSection("coordinates");
    234:       const Obj<typename Mesh::real_section_type>& parallelCoordinates = newMesh->getRealSection("coordinates");

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

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

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

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

    276:       PETSc::Log::Event("DistributeMesh").begin();
    277:       // Create the cell partition
    278:       Partitioner::createPartitionV(mesh, cellPartition, height);
    279:       if (mesh->debug()) {
    280:         PetscViewer    viewer;

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

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

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

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

    362:         for(std::set<std::string>::const_iterator n_iter = names->begin(); n_iter != names->end(); ++n_iter) {
    363:           const Obj<typename Mesh::int_section_type>& section    = mesh->getIntSection(*n_iter);
    364:           const Obj<typename Mesh::int_section_type>& newSection = newMesh->getIntSection(*n_iter);

    366:           // We assume all integer sections are complete sections
    367:           newSection->setChart(newMesh->getSieve()->getChart());
    368:           distributeSection(section, partition, renumbering, sendMeshOverlap, recvMeshOverlap, newSection);
    369:         }
    370:         PETSc::Log::Event("DistributeIntSec").end();
    371:       }
    372:       if (mesh->getArrowSections()->size() > 1) {
    373:         throw ALE::Exception("Need to distribute more arrow sections");
    374:       }
    375:       // Distribute labels
    376:       PETSc::Log::Event("DistributeLabels").begin();
    377:       const typename Mesh::labels_type& labels = mesh->getLabels();

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

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

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

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

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

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

    454:         ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, cones, overlapCones);
    455:         if (l->debug()) {overlapCones->view("Overlap Label Values");}
    456:         // Inserts cones into newL (must renumber here)
    457:         //ALE::Pullback::InsertionBinaryFusion::fuse(overlapCones, recvOverlap, renumbering, newSieve);
    458:         {
    459:           typedef typename cones_type::point_type overlap_point_type;
    460:           const typename RecvOverlap::capSequence::iterator rBegin = recvOverlap->capBegin();
    461:           const typename RecvOverlap::capSequence::iterator rEnd   = recvOverlap->capEnd();

    463:           for(typename RecvOverlap::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
    464:             const int                                             rank    = *r_iter;
    465:             const typename RecvOverlap::supportSequence::iterator pBegin  = recvOverlap->supportBegin(*r_iter);
    466:             const typename RecvOverlap::supportSequence::iterator pEnd    = recvOverlap->supportEnd(*r_iter);

    468:             for(typename RecvOverlap::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
    469:               const typename RecvOverlap::target_type& localPoint   = *p_iter;
    470:               const typename RecvOverlap::target_type& remotePoint  = p_iter.color();
    471:               const overlap_point_type                 overlapPoint = overlap_point_type(rank, remotePoint);
    472:               const int                                size         = overlapCones->getFiberDimension(overlapPoint);
    473:               const typename cones_type::value_type   *values       = overlapCones->restrictPoint(overlapPoint);

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

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

    522:       newMesh->setSieve(newSieve);
    523:       const Obj<partition_type> partition = unifyMesh(mesh, newMesh, renumbering, sendMeshOverlap, recvMeshOverlap);
    524:       // Unify coordinates
    525:       const Obj<typename Mesh::real_section_type>& coordinates    = mesh->getRealSection("coordinates");
    526:       const Obj<typename Mesh::real_section_type>& newCoordinates = newMesh->getRealSection("coordinates");

    528:       newMesh->setupCoordinates(newCoordinates);
    529:       distributeSection(coordinates, partition, renumbering, sendMeshOverlap, recvMeshOverlap, newCoordinates);
    530:       // Unify labels
    531:       const typename Mesh::labels_type& labels = mesh->getLabels();

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

    538:         //completeCones(label, newLabel, renumbering, sendMeshOverlap, recvMeshOverlap);
    539:         {
    540:           typedef ALE::UniformSection<point_type, int> cones_type;
    541:           typedef ALE::LabelSection<typename Mesh::sieve_type,typename Mesh::label_type> cones_wrapper_type;
    542:           Obj<cones_wrapper_type> cones        = new cones_wrapper_type(mesh->getSieve(), label);
    543:           Obj<cones_type>         overlapCones = new cones_type(label->comm(), label->debug());

    545:           ALE::Pullback::SimpleCopy::copy(sendMeshOverlap, recvMeshOverlap, cones, overlapCones);
    546:           if (label->debug()) {overlapCones->view("Overlap Label Values");}
    547:           // Inserts cones into parallelMesh (must renumber here)
    548:           //ALE::Pullback::InsertionBinaryFusion::fuse(overlapCones, recvMeshOverlap, renumbering, newSieve);
    549:           {
    550:             const Obj<typename mesh_recv_overlap_type::baseSequence> rPoints = recvMeshOverlap->base();

    552:             for(typename mesh_recv_overlap_type::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rPoints->end(); ++p_iter) {
    553:               const Obj<typename mesh_recv_overlap_type::coneSequence>& points      = recvMeshOverlap->cone(*p_iter);
    554:               const typename mesh_recv_overlap_type::target_type&       localPoint  = *p_iter;
    555:               const typename cones_type::point_type&                    remotePoint = points->begin().color();
    556:               const int                                                 size        = overlapCones->getFiberDimension(remotePoint);
    557:               const typename cones_type::value_type                    *values      = overlapCones->restrictPoint(remotePoint);

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

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

    648:           part_type *assignment = scatterBundle<Partitioner>(bundle, dim, bundleNew, sendOverlap, recvOverlap, height);
    649:           if (!subBundle.isNull() && !subBundleNew.isNull()) {
    650:             part_type *subAssignment = GenPartitioner::subordinatePartition(bundle, 1, subBundle, assignment);
    651:             const Obj<sieve_type>& sieve      = subBundle->getSieve();
    652:             const Obj<sieve_type>& sieveNew   = new typename ALE::Mesh<PetscInt,PetscScalar>::sieve_type(subBundle->comm(), subBundle->debug());
    653:             const int              numCells   = subBundle->heightStratum(height)->size();

    655:             subBundleNew->setSieve(sieveNew);
    656:             sieveCompletion::scatterSieve(subBundle, sieve, dim, sieveNew, sendOverlap, recvOverlap, height, numCells, subAssignment);
    657:             subBundleNew->stratify();
    658:             if (subAssignment != NULL) delete [] subAssignment;
    659:           }
    660:           if (assignment != NULL) delete [] assignment;
    661: #else
    662:           throw ALE::Exception("Chaco is not installed. Reconfigure with the flag --download-chaco");
    663: #endif
    664:         } else if (partitioner == "parmetis") {
    665: #ifdef PETSC_HAVE_PARMETIS
    666:           typedef typename ALE::New::ParMetis::Partitioner<bundle_type> Partitioner;
    667:           typedef typename ALE::New::Partitioner<bundle_type>           GenPartitioner;
    668:           typedef typename Partitioner::part_type                       part_type;

    670:           part_type *assignment = scatterBundle<Partitioner>(bundle, dim, bundleNew, sendOverlap, recvOverlap, height);
    671:           if (!subBundle.isNull() && !subBundleNew.isNull()) {
    672:             part_type *subAssignment = GenPartitioner::subordinatePartition(bundle, 1, subBundle, assignment);
    673:             const Obj<sieve_type>& sieve      = subBundle->getSieve();
    674:             const Obj<sieve_type>& sieveNew   = new typename ALE::Mesh<PetscInt,PetscScalar>::sieve_type(subBundle->comm(), subBundle->debug());
    675:             const int              numCells   = subBundle->heightStratum(height)->size();

    677:             subBundleNew->setSieve(sieveNew);
    678:             sieveCompletion::scatterSieve(subBundle, sieve, dim, sieveNew, sendOverlap, recvOverlap, height, numCells, subAssignment);
    679:             subBundleNew->stratify();
    680:             if (subAssignment != NULL) delete [] subAssignment;
    681:           }
    682:           if (assignment != NULL) delete [] assignment;
    683: #else
    684:           throw ALE::Exception("ParMetis is not installed. Reconfigure with the flag --download-parmetis");
    685: #endif
    686:         } else {
    687:           throw ALE::Exception("Unknown partitioner");
    688:         }
    689:       } else if (height == 1) {
    690:         if (partitioner == "zoltan") {
    691: #ifdef PETSC_HAVE_ZOLTAN
    692:           typedef typename ALE::New::Zoltan::Partitioner<bundle_type> Partitioner;
    693:           typedef typename Partitioner::part_type                     part_type;

    695:           part_type *assignment = scatterBundle<Partitioner>(bundle, dim, bundleNew, sendOverlap, recvOverlap, height);
    696:           if (assignment != NULL) delete [] assignment;
    697: #else
    698:           throw ALE::Exception("Zoltan is not installed. Reconfigure with the flag --download-zoltan");
    699: #endif
    700:         } else if (partitioner == "parmetis") {
    701: #ifdef PETSC_HAVE_PARMETIS
    702:           typedef typename ALE::New::ParMetis::Partitioner<bundle_type> Partitioner;
    703:           typedef typename Partitioner::part_type                       part_type;

    705:           part_type *assignment = scatterBundle<Partitioner>(bundle, dim, bundleNew, sendOverlap, recvOverlap, height);
    706:           if (assignment != NULL) delete [] assignment;
    707: #else
    708:           throw ALE::Exception("ParMetis is not installed. Reconfigure with the flag --download-parmetis");
    709: #endif
    710:         } else {
    711:           throw ALE::Exception("Unknown partitioner");
    712:         }
    713:       } else {
    714:         throw ALE::Exception("Invalid partition height");
    715:       }
    716:     }
    717:     template<typename Partitioner>
    718:     static typename Partitioner::part_type *scatterBundle(const Obj<bundle_type>& bundle, const int dim, const Obj<bundle_type>& bundleNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const int height = 0) {
    719:       typename Partitioner::part_type *assignment = createAssignment<Partitioner>(bundle, dim, sendOverlap, recvOverlap, height);
    720:       const Obj<sieve_type>&           sieve      = bundle->getSieve();
    721:       const Obj<sieve_type>&           sieveNew   = bundleNew->getSieve();
    722:       const int                        numPoints  = bundle->heightStratum(height)->size();

    724:       sieveCompletion::scatterSieve(bundle, sieve, dim, sieveNew, sendOverlap, recvOverlap, height, numPoints, assignment);
    725:       bundleNew->stratify();
    726:       return assignment;
    727:     }
    730:     static Obj<ALE::Mesh<PetscInt,PetscScalar> > distributeMesh(const Obj<ALE::Mesh<PetscInt,PetscScalar> >& serialMesh, const int height = 0, const std::string& partitioner = "chaco") {
    731:       typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
    732:       MPI_Comm                         comm          = serialMesh->comm();
    733:       const int                        dim           = serialMesh->getDimension();
    734:       Obj<FlexMesh>                    parallelMesh  = new FlexMesh(comm, dim, serialMesh->debug());
    735:       const Obj<FlexMesh::sieve_type>& parallelSieve = new FlexMesh::sieve_type(comm, serialMesh->debug());

    737:       ALE_LOG_EVENT_BEGIN;
    738:       parallelMesh->setSieve(parallelSieve);
    739:       if (serialMesh->debug()) {serialMesh->view("Serial mesh");}

    741:       // Distribute cones
    742:       Obj<send_overlap_type> sendOverlap = new send_overlap_type(comm, serialMesh->debug());
    743:       Obj<recv_overlap_type> recvOverlap = new recv_overlap_type(comm, serialMesh->debug());
    744:       scatterBundle(serialMesh, dim, parallelMesh, sendOverlap, recvOverlap, partitioner, height);
    745:       parallelMesh->setDistSendOverlap(sendOverlap);
    746:       parallelMesh->setDistRecvOverlap(recvOverlap);

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

    751:       for(typename bundle_type::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
    752:         if (parallelMesh->hasLabel(l_iter->first)) continue;
    753:         const Obj<typename bundle_type::label_type>& serialLabel   = l_iter->second;
    754:         const Obj<typename bundle_type::label_type>& parallelLabel = parallelMesh->createLabel(l_iter->first);
    755:         // Create local label
    756: #define NEW_LABEL
    757: #ifdef NEW_LABEL
    758:         parallelLabel->add(serialLabel, parallelSieve);
    759: #else
    760:         const Obj<typename bundle_type::label_type::traits::baseSequence>& base = serialLabel->base();

    762:         for(typename bundle_type::label_type::traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
    763:           if (parallelSieve->capContains(*b_iter) || parallelSieve->baseContains(*b_iter)) {
    764:             parallelLabel->addArrow(*serialLabel->cone(*b_iter)->begin(), *b_iter);
    765:           }
    766:         }
    767: #endif
    768:         // Get remote labels
    769:         sieveCompletion::scatterCones(serialLabel, parallelLabel, sendOverlap, recvOverlap);
    770:       }

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

    775:       for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
    776:         parallelMesh->setRealSection(*name, distributeSection(serialMesh->getRealSection(*name), parallelMesh, sendOverlap, recvOverlap));
    777:       }
    778:       sections = serialMesh->getIntSections();
    779:       for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
    780:         parallelMesh->setIntSection(*name, distributeSection(serialMesh->getIntSection(*name), parallelMesh, sendOverlap, recvOverlap));
    781:       }
    782:       sections = serialMesh->getArrowSections();

    784:       for(std::set<std::string>::iterator name = sections->begin(); name != sections->end(); ++name) {
    785:         parallelMesh->setArrowSection(*name, distributeArrowSection(serialMesh->getArrowSection(*name), serialMesh, parallelMesh, sendOverlap, recvOverlap));
    786:       }
    787:       if (parallelMesh->debug()) {parallelMesh->view("Parallel Mesh");}
    788:       ALE_LOG_EVENT_END;
    789:       return parallelMesh;
    790:     }
    793:     template<typename Section>
    794:     static void updateSectionLocal(const Obj<Section>& oldSection, const Obj<bundle_type>& newBundle, const Obj<Section>& newSection) {
    795:       const Obj<typename bundle_type::sieve_type>&    newSieve = newBundle->getSieve();
    796:       const typename Section::atlas_type::chart_type& oldChart = oldSection->getChart();

    798:       for(typename Section::atlas_type::chart_type::const_iterator c_iter = oldChart.begin(); c_iter != oldChart.end(); ++c_iter) {
    799:         if (newSieve->capContains(*c_iter) || newSieve->baseContains(*c_iter)) {
    800:           newSection->setFiberDimension(*c_iter, oldSection->getFiberDimension(*c_iter));
    801:         }
    802:       }
    803:       newBundle->allocate(newSection);
    804:       const typename Section::atlas_type::chart_type& newChart = newSection->getChart();

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

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

    820:         for(typename recv_overlap_type::traits::coneSequence::iterator p_iter = recvPatches->begin(); p_iter != end; ++p_iter) {
    821:           newSection->addPoint(*r_iter, recvSection->getSection(*p_iter)->getFiberDimension(*r_iter));
    822:         }
    823:       }
    824:       newBundle->reallocate(newSection);
    825:       for(typename recv_overlap_type::traits::baseSequence::iterator r_iter = recvPoints->begin(); r_iter != recvPoints->end(); ++r_iter) {
    826:         const Obj<typename recv_overlap_type::traits::coneSequence>&     recvPatches = recvOverlap->cone(*r_iter);
    827:         const typename recv_overlap_type::traits::coneSequence::iterator end         = recvPatches->end();

    829:         for(typename recv_overlap_type::traits::coneSequence::iterator p_iter = recvPatches->begin(); p_iter != end; ++p_iter) {
    830:           if (recvSection->getSection(*p_iter)->getFiberDimension(*r_iter)) {
    831:             newSection->updatePointAll(*r_iter, recvSection->getSection(*p_iter)->restrictPoint(*r_iter));
    832:           }
    833:         }
    834:       }
    835:     }
    838:     template<typename Section>
    839:     static Obj<Section> distributeSection(const Obj<Section>& serialSection, const Obj<bundle_type>& parallelBundle, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
    840:       if (serialSection->debug()) {
    841:         serialSection->view("Serial Section");
    842:       }
    843:       typedef typename alloc_type::template rebind<typename Section::value_type>::other value_alloc_type;
    844:       typedef ALE::Field<send_overlap_type, int, ALE::Section<point_type, typename Section::value_type, value_alloc_type> > send_section_type;
    845:       typedef ALE::Field<recv_overlap_type, int, ALE::Section<point_type, typename Section::value_type, value_alloc_type> > recv_section_type;
    846:       typedef ALE::New::SizeSection<Section> SectionSizer;
    847:       Obj<Section>                 parallelSection = new Section(serialSection->comm(), serialSection->debug());
    848:       const Obj<send_section_type> sendSection     = new send_section_type(serialSection->comm(), serialSection->debug());
    849:       const Obj<recv_section_type> recvSection     = new recv_section_type(serialSection->comm(), sendSection->getTag(), serialSection->debug());
    850:       const Obj<SectionSizer>      sizer           = new SectionSizer(serialSection);

    852:       updateSectionLocal(serialSection, parallelBundle, parallelSection);
    853:       sectionCompletion::completeSection(sendOverlap, recvOverlap, sizer, serialSection, sendSection, recvSection);
    854:       updateSectionRemote(recvOverlap, recvSection, parallelBundle, parallelSection);
    855:       if (parallelSection->debug()) {
    856:         parallelSection->view("Parallel Section");
    857:       }
    858:       return parallelSection;
    859:     }
    862:     template<typename Section>
    863:     static void updateArrowSectionLocal(const Obj<Section>& oldSection, const Obj<bundle_type>& newBundle, const Obj<Section>& newSection) {
    864:       const Obj<typename bundle_type::sieve_type>&    newSieve = newBundle->getSieve();
    865:       const typename Section::atlas_type::chart_type& oldChart = oldSection->getChart();

    867:       for(typename Section::atlas_type::chart_type::const_iterator c_iter = oldChart.begin(); c_iter != oldChart.end(); ++c_iter) {
    868:         // Dmitry should provide a Sieve::contains(MinimalArrow) method
    869:         if (newSieve->capContains(c_iter->source) && newSieve->baseContains(c_iter->target)) {
    870:           newSection->setFiberDimension(*c_iter, oldSection->getFiberDimension(*c_iter));
    871:         }
    872:       }
    873:       //newBundle->allocate(newSection);
    874:       const typename Section::atlas_type::chart_type& newChart = newSection->getChart();

    876:       for(typename Section::atlas_type::chart_type::const_iterator c_iter = newChart.begin(); c_iter != newChart.end(); ++c_iter) {
    877:         newSection->updatePointAll(*c_iter, oldSection->restrictPoint(*c_iter));
    878:       }
    879:     }
    882:     template<typename RecvSection, typename Section>
    883:     static void updateArrowSectionRemote(const Obj<recv_overlap_type>& recvOverlap, const Obj<RecvSection>& recvSection, const Obj<bundle_type>& newBundle, const Obj<Section>& newSection) {
    884:       Obj<typename recv_overlap_type::traits::baseSequence> recvPoints = recvOverlap->base();

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

    890:         for(typename bundle_type::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != end; ++c_iter) {
    891:           newSection->setFiberDimension(typename Section::point_type(*c_iter, *r_iter), 1);
    892:         }
    893:       }
    894:       //newBundle->reallocate(newSection);
    895:       for(typename recv_overlap_type::traits::baseSequence::iterator r_iter = recvPoints->begin(); r_iter != recvPoints->end(); ++r_iter) {
    896:         const Obj<typename recv_overlap_type::traits::coneSequence>&     recvPatches = recvOverlap->cone(*r_iter);
    897:         const typename recv_overlap_type::traits::coneSequence::iterator recvEnd     = recvPatches->end();

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

    902:           if (section->getFiberDimension(*r_iter)) {
    903:             const Obj<typename bundle_type::sieve_type::traits::coneSequence>&     cone    = newBundle->getSieve()->cone(*r_iter);
    904:             const typename bundle_type::sieve_type::traits::coneSequence::iterator end     = cone->end();
    905:             const typename RecvSection::value_type                                *values  = section->restrictPoint(*r_iter);
    906:             int                                                                    c       = -1;

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

    933:       updateArrowSectionLocal(serialSection, parallelBundle, parallelSection);
    934:       sectionCompletion::completeSection(sendOverlap, recvOverlap, sizer, filler, sendSection, recvSection);
    935:       updateArrowSectionRemote(recvOverlap, recvSection, parallelBundle, parallelSection);
    936:       if (parallelSection->debug()) {
    937:         parallelSection->view("Parallel ArrowSection");
    938:       }
    939:       return parallelSection;
    940:     }
    941:     static void unifyBundle(const Obj<bundle_type>& bundle, const int dim, const Obj<bundle_type>& bundleNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
    942:       typedef int part_type;
    943:       const Obj<sieve_type>& sieve    = bundle->getSieve();
    944:       const Obj<sieve_type>& sieveNew = bundleNew->getSieve();
    945:       const int              rank     = bundle->commRank();
    946:       const int              debug    = bundle->debug();

    948:       // 1) Form partition point overlap a priori
    949:       if (rank == 0) {
    950:         for(int p = 1; p < sieve->commSize(); p++) {
    951:           // The arrow is from remote partition point 0 on rank p to local partition point 0
    952:           recvOverlap->addCone(p, 0, 0);
    953:         }
    954:       } else {
    955:         // The arrow is from local partition point 0 to remote partition point 0 on rank 0
    956:         sendOverlap->addCone(0, 0, 0);
    957:       }
    958:       if (debug) {
    959:         sendOverlap->view("Send overlap for partition");
    960:         recvOverlap->view("Receive overlap for partition");
    961:       }
    962:       // 2) Partition the mesh
    963:       int        numCells   = bundle->heightStratum(0)->size();
    964:       part_type *assignment = new part_type[numCells];

    966:       for(int c = 0; c < numCells; ++c) {
    967:         assignment[c] = 0;
    968:       }
    969:       // 3) Scatter the sieve
    970:       sieveCompletion::scatterSieve(bundle, sieve, dim, sieveNew, sendOverlap, recvOverlap, 0, numCells, assignment);
    971:       bundleNew->stratify();
    972:       // 4) Cleanup
    973:       if (assignment != NULL) delete [] assignment;
    974:     }
    977:     static Obj<ALE::Mesh<PetscInt,PetscScalar> > unifyMesh(const Obj<ALE::Mesh<PetscInt,PetscScalar> >& parallelMesh) {
    978:       typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
    979:       const int                        dim         = parallelMesh->getDimension();
    980:       Obj<FlexMesh>                    serialMesh  = new FlexMesh(parallelMesh->comm(), dim, parallelMesh->debug());
    981:       const Obj<FlexMesh::sieve_type>& serialSieve = new FlexMesh::sieve_type(parallelMesh->comm(), parallelMesh->debug());

    983:       ALE_LOG_EVENT_BEGIN;
    984:       serialMesh->setSieve(serialSieve);
    985:       if (parallelMesh->debug()) {
    986:         parallelMesh->view("Parallel topology");
    987:       }

    989:       // Unify cones
    990:       Obj<send_overlap_type> sendOverlap = new send_overlap_type(serialMesh->comm(), serialMesh->debug());
    991:       Obj<recv_overlap_type> recvOverlap = new recv_overlap_type(serialMesh->comm(), serialMesh->debug());
    992:       unifyBundle(parallelMesh, dim, serialMesh, sendOverlap, recvOverlap);
    993:       serialMesh->setDistSendOverlap(sendOverlap);
    994:       serialMesh->setDistRecvOverlap(recvOverlap);

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

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

    1004:         sieveCompletion::scatterCones(parallelLabel, serialLabel, sendOverlap, recvOverlap);
    1005:       }

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

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

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

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

    1044:           for(int p = 0; p < size; p++) {
    1045:             if (origSendOverlap->support(points[p])->size() == 0) {
    1046:               sendOverlap->addArrow(points[p], rank, points[p]);
    1047:             }
    1048:           }
    1049:         }
    1050:       }
    1051:       for(typename RecvSection::sheaf_type::const_iterator p_iter = recvRanks.begin(); p_iter != recvRanks.end(); ++p_iter) {
    1052:         const typename RecvSection::patch_type&               rank    = p_iter->first;
    1053:         const Obj<typename RecvSection::section_type>&        section = p_iter->second;
    1054:         const typename RecvSection::section_type::chart_type& chart   = section->getChart();

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

    1060:           for(int p = 0; p < size; p++) {
    1061:             if (origRecvOverlap->support(rank, points[p])->size() == 0) {
    1062:               recvOverlap->addArrow(rank, points[p], points[p]);
    1063:             }
    1064:           }
    1065:         }
    1066:       }
    1067:     }
    1070:     template<typename RecvOverlap, typename RecvSection>
    1071:     static void updateSieve(const Obj<RecvOverlap>& recvOverlap, const Obj<RecvSection>& recvSection, const Obj<sieve_type>& sieve) {
    1072: #if 1
    1073:       Obj<typename RecvOverlap::traits::baseSequence> recvPoints = recvOverlap->base();

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

    1079:         for(typename RecvOverlap::traits::coneSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
    1080:           const typename RecvOverlap::target_type&       remotePoint = r_iter.color();
    1081:           const int                                      rank        = *r_iter;
    1082:           const Obj<typename RecvSection::section_type>& section     = recvSection->getSection(rank);
    1083:           const typename RecvSection::value_type        *points      = section->restrictPoint(remotePoint);
    1084:           const int                                      size        = section->getFiberDimension(remotePoint);
    1085:           int                                            c           = 0;

    1087:           ///std::cout << "["<<recvSection->commRank()<<"]: Receiving " << size << " points from rank " << rank << std::endl;
    1088:           for(int p = 0; p < size; p++) {
    1089:             // rank -- remote point --> local point
    1090:             if (recvOverlap->support(rank, points[p])->size()) {
    1091:               sieve->addArrow(*recvOverlap->support(rank, points[p])->begin(), localPoint, c);
    1092:               ///std::cout << "["<<recvSection->commRank()<<"]:   1Adding arrow " << *recvOverlap->support(rank, points[p])->begin() << "("<<points[p]<<") --> " << localPoint << std::endl;
    1093:             } else {
    1094:               sieve->addArrow(points[p], localPoint, c);
    1095:               ///std::cout << "["<<recvSection->commRank()<<"]:   2Adding arrow " << points[p] << " --> " << localPoint << std::endl;
    1096:             }
    1097:           }
    1098:         }
    1099:       }
    1100: #else
    1101:       const typename RecvSection::sheaf_type& ranks = recvSection->getPatches();

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

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

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

    1150:       bundle->constructOverlap();
    1151:       const Obj<send_overlap_type> sendOverlap = bundle->getSendOverlap();
    1152:       const Obj<recv_overlap_type> recvOverlap = bundle->getRecvOverlap();
    1153:       const Obj<send_section_type> sendSection = new send_section_type(section->comm(), section->debug());
    1154:       const Obj<recv_section_type> recvSection = new recv_section_type(section->comm(), sendSection->getTag(), section->debug());
    1155:       const Obj<SectionSizer>      sizer       = new SectionSizer(section);

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

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

    1165:         for(typename recv_overlap_type::traits::coneSequence::iterator p_iter = recvPatches->begin(); p_iter != end; ++p_iter) {
    1166:           if (recvSection->getSection(*p_iter)->getFiberDimension(p_iter.color())) {
    1167:             if (debug) {std::cout << "["<<section->commRank()<<"]Completed point " << *r_iter << std::endl;}
    1168:             section->updateAddPoint(*r_iter, recvSection->getSection(*p_iter)->restrictPoint(p_iter.color()));
    1169:           }
    1170:         }
    1171:       }
    1172:     }
    1173:   };
    1174: }
    1175: #endif