Commits

Matt Knepley committed 6aa6582

Lots of work on Sieve
- Nearly have preallocation done

Hg-commit: fa2bb7a04343cbd69c5997bc1b38860e994d407c

  • Participants
  • Parent commits eedeaf5

Comments (0)

Files changed (7)

src/dm/mesh/examples/tests/mesh1.cxx

   const Obj<ALE::Mesh::numbering_type>& order   = mesh->getGlobalOrder(name);
   const Obj<ALE::Mesh::section_type>&   section = mesh->getSection(name);
   Mat                                   A;
+  PetscInt                              numLocalRows, firstRow, lastRow;
+  PetscInt                             *dnz, *onz;
   PetscErrorCode                        ierr;
 
   PetscFunctionBegin;
   ierr = MatCreate(mesh->comm(), &A);CHKERRQ(ierr);
   ierr = MatSetSizes(A, order->getLocalSize(), order->getLocalSize(), order->getGlobalSize(), order->getGlobalSize());CHKERRQ(ierr);
   ierr = MatSetFromOptions(A);CHKERRQ(ierr);
-  ierr = preallocateMatrix(mesh.objPtr, section, order, A);CHKERRQ(ierr);
+  //ierr = preallocateMatrix(mesh.objPtr, section, order, A);CHKERRQ(ierr);
+  const Obj<ALE::Mesh::sieve_type>    adjGraph    = new ALE::Mesh::sieve_type(mesh->comm(), mesh->debug);
+  const Obj<ALE::Mesh::topology_type> adjTopology = new ALE::Mesh::topology_type(mesh->comm(), mesh->debug);
+  const ALE::Mesh::section_type::patch_type patch = 0;
+
+  adjTopology->setPatch(patch, adjGraph);
+  ierr = MatGetLocalSize(A, &numLocalRows, PETSC_NULL);CHKERRQ(ierr);
+  ierr = MatGetOwnershipRange(A, &firstRow, &lastRow);CHKERRQ(ierr);
+  ierr = PetscMalloc2(numLocalRows, PetscInt, &dnz, numLocalRows, PetscInt, &onz);CHKERRQ(ierr);
+  // Create local adjacency graph
+  //   In general, we need to FIAT inof that attaches dual basis vectors to sieve points
+  const ALE::Mesh::atlas_type::chart_type& chart = section->getAtlas()->getChart(patch);
+  const Obj<ALE::Mesh::sieve_type>&        sieve = mesh->getTopologyNew()->getPatch(patch);
+
+  for(ALE::Mesh::atlas_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
+    const ALE::Mesh::atlas_type::point_type& point = c_iter->first;
+
+    adjGraph->addCone(sieve->cone(sieve->support(point)), point);
+  }
+  // Distribute adjacency graph
+  const Obj<ALE::Mesh::numbering_type>&    vNumbering        = mesh->getNumbering(0);
+  const Obj<ALE::Mesh::send_overlap_type>& vertexSendOverlap = vNumbering->getSendOverlap();
+  const Obj<ALE::Mesh::recv_overlap_type>& vertexRecvOverlap = vNumbering->getRecvOverlap();
+  const Obj<ALE::Mesh::send_overlap_type>  nbrSendOverlap    = new ALE::Mesh::send_overlap_type(mesh->comm(), mesh->debug);
+  const Obj<ALE::Mesh::recv_overlap_type>  nbrRecvOverlap    = new ALE::Mesh::recv_overlap_type(mesh->comm(), mesh->debug);
+  const Obj<ALE::Mesh::send_section_type>  sendSection       = new ALE::Mesh::send_section_type(mesh->comm(), mesh->debug);
+  const Obj<ALE::Mesh::recv_section_type>  recvSection       = new ALE::Mesh::recv_section_type(mesh->comm(), sendSection->getTag(), mesh->debug);
+
+  ALE::New::Distribution<ALE::Mesh::topology_type>::coneCompletion(vertexSendOverlap, vertexRecvOverlap, adjTopology, sendSection, recvSection);
+  adjTopology->view("Adjacency topology");
+  vNumbering->view("Global vertex numbering");
+  order->view("Global vertex order");
+  // Distribute indices for new points
+  ALE::New::Distribution<ALE::Mesh::topology_type>::updateOverlap(sendSection, recvSection, nbrSendOverlap, nbrRecvOverlap);
+  order->setSendOverlap(nbrSendOverlap);
+  order->setRecvOverlap(nbrRecvOverlap);
+  nbrSendOverlap->view("Neighbor Send Overlap");
+  nbrRecvOverlap->view("Neighbor Receive Overlap");
+  order->complete(true);
+  order->view("Global vertex order after completion");
+  // Read out adjacency graph
+  const ALE::Obj<ALE::Mesh::sieve_type> graph = adjTopology->getPatch(patch);
+
+  ierr = PetscMemzero(dnz, numLocalRows * sizeof(PetscInt));CHKERRQ(ierr);
+  ierr = PetscMemzero(onz, numLocalRows * sizeof(PetscInt));CHKERRQ(ierr);
+  for(ALE::Mesh::atlas_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
+    const ALE::Mesh::atlas_type::point_type& point = c_iter->first;
+
+    if (order->isLocal(point)) {
+      const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& adj   = graph->cone(point);
+      const int                                               row   = order->getIndex(point);
+      const int                                               rSize = section->getAtlas()->getFiberDimension(patch, point);
+
+      for(ALE::Mesh::sieve_type::traits::coneSequence::iterator v_iter = adj->begin(); v_iter != adj->end(); ++v_iter) {
+        const ALE::Mesh::atlas_type::point_type& neighbor = *v_iter;
+        const int                                col      = order->getIndex(neighbor);
+        const int                                cSize    = section->getAtlas()->getFiberDimension(patch, neighbor);
+        
+        if (col >= firstRow && col < lastRow) {
+          for(int r = 0; r < rSize; ++r) {dnz[row - firstRow + r] += cSize;}
+        } else {
+          for(int r = 0; r < rSize; ++r) {onz[row - firstRow + r] += cSize;}
+        }
+      }
+    }
+  }
+  int rank = mesh->commRank();
+  for(int r = 0; r < numLocalRows; r++) {
+    std::cout << "["<<rank<<"]: dnz["<<r<<"]: " << dnz[r] << " onz["<<r<<"]: " << onz[r] << std::endl;
+  }
+  ierr = MatSeqAIJSetPreallocation(A, 0, dnz);CHKERRQ(ierr);
+  ierr = MatMPIAIJSetPreallocation(A, 0, dnz, 0, onz);CHKERRQ(ierr);
+  ierr = PetscFree2(dnz, onz);CHKERRQ(ierr);
+  ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR);CHKERRQ(ierr);
+  // Fill matrix
   const Obj<ALE::Mesh::topology_type::label_sequence>& elements = mesh->getTopologyNew()->heightStratum(0, 0);
   int          size       = options->dim*options->dim*9;
   PetscScalar *elementMat = new PetscScalar[size];

src/dm/mesh/mesh.c

 #define __FUNCT__ "preallocateMatrix"
 PetscErrorCode preallocateMatrix(ALE::Mesh *mesh, const ALE::Obj<ALE::Mesh::section_type>& field, const ALE::Obj<ALE::Mesh::numbering_type>& globalOrder, Mat A)
 {
-        ALE::Obj<ALE::Mesh::topology_type>  parallelTopology = new ALE::Mesh::topology_type(mesh->comm(), mesh->debug);
-  const ALE::Obj<ALE::Mesh::sieve_type>     serialGraph      = new ALE::Mesh::sieve_type(mesh->comm(), mesh->debug);
-  const ALE::Obj<ALE::Mesh::topology_type>  serialTopology   = new ALE::Mesh::topology_type(mesh->comm(), mesh->debug);
-  const ALE::Mesh::section_type::patch_type patch = 0;
-  const ALE::Obj<ALE::Mesh::sieve_type>&    sieve = mesh->getTopologyNew()->getPatch(patch);
+  const ALE::Obj<ALE::Mesh::sieve_type>     adjGraph    = new ALE::Mesh::sieve_type(mesh->comm(), mesh->debug);
+  const ALE::Obj<ALE::Mesh::topology_type>  adjTopology = new ALE::Mesh::topology_type(mesh->comm(), mesh->debug);
+  const ALE::Mesh::section_type::patch_type patch       = 0;
+  const ALE::Obj<ALE::Mesh::sieve_type>&    sieve       = mesh->getTopologyNew()->getPatch(patch);
   PetscInt       numLocalRows, firstRow, lastRow;
   PetscInt      *dnz, *onz;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
-  serialTopology->setPatch(patch, serialGraph);
+  adjTopology->setPatch(patch, adjGraph);
   ierr = MatGetLocalSize(A, &numLocalRows, PETSC_NULL);CHKERRQ(ierr);
   ierr = MatGetOwnershipRange(A, &firstRow, &lastRow);CHKERRQ(ierr);
   ierr = PetscMalloc2(numLocalRows, PetscInt, &dnz, numLocalRows, PetscInt, &onz);CHKERRQ(ierr);
   for(ALE::Mesh::atlas_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
     const ALE::Mesh::atlas_type::point_type& point = c_iter->first;
 
-    serialGraph->addCone(sieve->cone(sieve->support(point)), point);
+    adjGraph->addCone(sieve->cone(sieve->support(point)), point);
   }
   /* Distribute adjacency graph */
-  const ALE::Obj<ALE::Mesh::send_overlap_type>&              vertexSendOverlap = mesh->getVertexSendOverlap();
-  const ALE::Obj<ALE::Mesh::recv_overlap_type>&              vertexRecvOverlap = mesh->getVertexRecvOverlap();
+#if 1
+  const ALE::Obj<ALE::Mesh::send_overlap_type>& vertexSendOverlap = mesh->getVertexSendOverlap();
+  const ALE::Obj<ALE::Mesh::recv_overlap_type>& vertexRecvOverlap = mesh->getVertexRecvOverlap();
+
+  //ALE::New::Distribution<ALE::Mesh::topology_type>::coneCompletion(vertexSendOverlap, vertexRecvOverlap, adjTopology);
+#else
+  const ALE::Obj<ALE::Mesh::send_overlap_type>&             vertexSendOverlap = mesh->getVertexSendOverlap();
+  const ALE::Obj<ALE::Mesh::recv_overlap_type>&             vertexRecvOverlap = mesh->getVertexRecvOverlap();
   const ALE::Obj<ALE::Mesh::numbering_type>&                vNumbering  = mesh->getLocalNumbering(0);
   const ALE::Obj<ALE::Mesh::topology_type::label_sequence>& vertices    = field->getAtlas()->getTopology()->depthStratum(patch, 0);
   int                                                       numVertices = vertices->size();
     }
   }
   ALE::New::Distribution<ALE::Mesh::topology_type>::scatterTopology(serialTopology, numVertices, assignment, parallelTopology);
+#endif
   /* Read out adjacency graph */
-  const ALE::Obj<ALE::Mesh::sieve_type> parallelGraph = parallelTopology->getPatch(patch);
+  const ALE::Obj<ALE::Mesh::sieve_type> graph = adjTopology->getPatch(patch);
 
   ierr = PetscMemzero(dnz, numLocalRows * sizeof(PetscInt));CHKERRQ(ierr);
   ierr = PetscMemzero(onz, numLocalRows * sizeof(PetscInt));CHKERRQ(ierr);
     const ALE::Mesh::atlas_type::point_type& point = c_iter->first;
 
     if (globalOrder->isLocal(point)) {
-      const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& adj   = parallelGraph->cone(point);
+      const Obj<ALE::Mesh::sieve_type::traits::coneSequence>& adj   = graph->cone(point);
       const int                                               row   = globalOrder->getIndex(point);
       const int                                               rSize = c_iter->second.index;
 
       }
     }
   }
+  int rank = mesh->commRank();
   for(int r = 0; r < numLocalRows; r++) {
     std::cout << "["<<rank<<"]: dnz["<<r<<"]: " << dnz[r] << " onz["<<r<<"]: " << onz[r] << std::endl;
   }

src/dm/mesh/sieve/CoSieve.hh

       const Obj<topology_type>& getTopology() const {return this->_topology;};
       void setTopology(const Obj<topology_type>& topology) {this->_topology = topology;};
       void copy(const Obj<Atlas>& atlas) {
-        const typename topology_type::sheaf_type& patches = atlas->getTopology()->getPatches();
+        const typename topology_type::sheaf_type& sheaf = atlas->getTopology()->getPatches();
 
-        for(typename topology_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
-          for(typename chart_type::const_iterator c_iter = p_iter->second.begin(); c_iter != p_iter->second.end(); ++c_iter) {
-            this->setFiberDimension(p_iter->first, c_iter->first, c_iter->second);
+        for(typename topology_type::sheaf_type::const_iterator s_iter = sheaf.begin(); s_iter != sheaf.end(); ++s_iter) {
+          const chart_type& chart = atlas->getChart(s_iter->first);
+
+          for(typename chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
+            this->setFiberDimension(s_iter->first, c_iter->first, c_iter->second.index);
           }
         }
-        this->orderPatches();
       };
       void copyByDepth(const Obj<Atlas>& atlas) {
         this->copyByDepth(atlas, atlas->getTopology());
       typedef std::map<patch_type, value_type *> values_type;
     protected:
       Obj<atlas_type> _atlas;
+      Obj<atlas_type> _atlasNew;
       values_type     _arrays;
+      typename atlas_type::indices_type _newPoints;
     public:
       Section(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
-        this->_atlas = new atlas_type(comm, debug);
+        this->_atlas    = new atlas_type(comm, debug);
+        this->_atlasNew = NULL;
       };
-      Section(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas) {};
+      Section(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas), _atlasNew(NULL) {};
       virtual ~Section() {
         for(typename values_type::iterator a_iter = this->_arrays.begin(); a_iter != this->_arrays.end(); ++a_iter) {
           delete [] a_iter->second;
       void updatePoint(const patch_type& patch, const point_type& p, const value_type v[]) {
         this->checkPatch(patch);
         const index_type& ind = this->_atlas->getIndex(patch, p);
-        value_type       *a   = &(this->_arrays[patch][ind.first]);
+        value_type       *a   = &(this->_arrays[patch][ind.prefix]);
 
         // Using the index structure explicitly
-        for(int i = 0; i < ind.second; ++i) {
+        for(int i = 0; i < ind.index; ++i) {
           a[i] = v[i];
         }
       };
           }
         }
       };
+    public: // Resizing
+      void addPoint(const patch_type& patch, const point_type& point, const int size) {
+        const typename atlas_type::chart_type& chart = this->_atlas->getChart(patch);
+
+        if (chart.find(point) == chart.end()) {
+          if (this->_atlasNew.isNull()) {
+            this->_atlasNew = new atlas_type(this->_atlas->getTopology());
+            this->_atlasNew->copy(this->_atlas);
+          }
+          this->_atlasNew->setFiberDimension(patch, point, size);
+        }
+      };
+      void reallocate() {
+        if (this->_atlasNew.isNull()) return;
+        this->_atlasNew->orderPatches();
+        const typename atlas_type::topology_type::sheaf_type& patches = this->_atlas->getTopology()->getPatches();
+
+        for(typename atlas_type::topology_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
+          const typename atlas_type::patch_type& patch    = p_iter->first;
+          const typename atlas_type::chart_type& chart    = this->_atlas->getChart(patch);
+          const value_type                      *array    = this->_arrays[patch];
+          value_type                            *newArray = new value_type[this->_atlasNew->size(patch)];
+
+          for(typename atlas_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
+            const int& size      = c_iter->second.index;
+            const int& offset    = c_iter->second.prefix;
+            const int& newOffset = this->_atlasNew->getIndex(patch, c_iter->first).prefix;
+
+            for(int i = 0; i < size; ++i) {
+              newArray[newOffset+i] = array[offset+i];
+            }
+          }
+          delete [] this->_arrays[patch];
+          this->_arrays[patch] = newArray;
+        }
+        this->_atlas    = this->_atlasNew;
+        this->_atlasNew = NULL;
+      };
     public:
       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
         ostringstream txt;

src/dm/mesh/sieve/Completion.hh

       typedef typename ALE::New::ConeSizeSection<topology_type, sieve_type>               cone_size_section;
       typedef typename ALE::New::ConeSection<topology_type, sieve_type>                   cone_section;
     public:
+      // Creates a DiscreteTopology with the overlap information
       static Obj<topology_type> createSendTopology(const Obj<send_overlap_type>& sendOverlap) {
         const Obj<send_overlap_type::traits::baseSequence> ranks = sendOverlap->base();
         Obj<topology_type> topology = new topology_type(sendOverlap->comm(), sendOverlap->debug);
         delete [] this->_offsets;
       };
     public: // Accessors
+      const Obj<topology_type>& getTopology() {return this->_atlas->getTopology();};
       std::string getLabel() const {return this->_label;};
       void        setLabel(const std::string& label) {this->_label = label;};
       int         getValue() const {return this->_value;};
       void        setValue(const int value) {this->_value = value;};
+      const Obj<send_overlap_type>& getSendOverlap() {return this->_sendOverlap;};
+      void setSendOverlap(const Obj<send_overlap_type>& overlap) {this->_sendOverlap = overlap;};
+      const Obj<recv_overlap_type>& getRecvOverlap() {return this->_recvOverlap;};
+      void setRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_recvOverlap = overlap;};
+      //const Obj<send_section_type>& getSendSection() {return this->_sendSection;};
+      //void setSendSection(const Obj<send_section_type>& section) {this->_sendSection = section;};
+      //const Obj<recv_section_type>& getRecvSection() {return this->_recvSection;};
+      //void setRecvSection(const Obj<recv_section_type>& section) {this->_recvSection = section;};
+      void copyCommunication(const Obj<NewNumbering<atlas_type> >& numbering) {
+        this->setSendOverlap(numbering->getSendOverlap());
+        this->setRecvOverlap(numbering->getRecvOverlap());
+        //this->setSendSection(numbering->getSendSection());
+        //this->setRecvSection(numbering->getRecvSection());
+      }
     public: // Sizes
       int        getLocalSize() const {return this->_localSize;};
+      void       setLocalSize(const int size) {this->_localSize = size;};
       int        getGlobalSize() const {return this->_offsets[this->commSize()];};
       const int *getGlobalOffsets() const {return this->_offsets;};
     public: // Indices
       int getIndex(const point_type& point) {if (this->restrictPoint(0, point)[0] >= 0) return this->restrictPoint(0, point)[0]; else return -(this->restrictPoint(0, point)[0]+1);};
+      void setIndex(const point_type& point, const int index) {this->updatePoint(0, point, &index);};
       point_type getPoint(const int& index) {return this->_invOrder[index];};
       bool isLocal(const point_type& point) {return this->restrictPoint(0, point)[0] >= 0;};
       bool isRemote(const point_type& point) {return this->restrictPoint(0, point)[0] < 0;};
+      bool hasPoint(const point_type& point) {
+        const typename atlas_type::chart_type& chart = this->getChart(0);
+        return (chart->find(point) != chart->end());
+      };
     public:
       void constructOverlap() {
         const Obj<typename topology_type::label_sequence>& points = this->getAtlas()->getTopology()->getLabelStratum(0, this->_label, this->_value);
           }
         }
       };
-      void complete() {
+      void complete(bool allowDuplicates = false) {
         typedef typename Completion<topology_type, int>::atlas_type atlas_type;
         typedef typename ALE::New::OverlapValues<send_overlap_type, atlas_type, value_type> send_section_type;
         typedef typename ALE::New::OverlapValues<recv_overlap_type, atlas_type, value_type> recv_section_type;
         Obj<typename recv_overlap_type::traits::baseSequence> recvPoints = this->_recvOverlap->base();
 
         for(typename recv_overlap_type::traits::baseSequence::iterator r_iter = recvPoints->begin(); r_iter != recvPoints->end(); ++r_iter) {
+          this->addPoint(0, *r_iter, 1);
+        }
+        this->reallocate();
+        for(typename recv_overlap_type::traits::baseSequence::iterator r_iter = recvPoints->begin(); r_iter != recvPoints->end(); ++r_iter) {
           const Obj<typename recv_overlap_type::traits::coneSequence>& recvPatches = this->_recvOverlap->cone(*r_iter);
     
           for(typename recv_overlap_type::traits::coneSequence::iterator p_iter = recvPatches->begin(); p_iter != recvPatches->end(); ++p_iter) {
             const typename recv_section_type::value_type *values = recvSection->restrict(*p_iter, *r_iter);
 
             if (values[0] >= 0) {
-              if (this->isLocal(*r_iter)) {
+              if (this->isLocal(*r_iter) && !allowDuplicates) {
                 ostringstream msg;
-                msg << "Multiple indices for point " << *r_iter;
+                msg << "["<<this->commRank()<<"]Multiple indices for point " << *r_iter << " from " << *p_iter << " with index " << values[0];
+                throw ALE::Exception(msg.str().c_str());
+              }
+              if (this->_atlas->getFiberDimension(0, *r_iter) == 0) {
+                ostringstream msg;
+                msg << "["<<this->commRank()<<"]Unexpected point " << *r_iter << " from " << *p_iter << " with index " << values[0];
                 throw ALE::Exception(msg.str().c_str());
               }
               int val = -(values[0]+1);
       Obj<recv_section_type>    _recvSection;
       int                       _localSize;
       int                      *_offsets;
-    public:
-      Numbering(const Obj<topology_type>& topology, const std::string& label, int value) : ParallelObject(topology->comm(), topology->debug()), _topology(topology), _label(label), _value(value), _localSize(0) {
+    protected:
+      void __init() {
         this->_sendOverlap = new send_overlap_type(this->comm(), this->debug());
         this->_recvOverlap = new recv_overlap_type(this->comm(), this->debug());
         this->_sendSection = new send_section_type(this->comm(), this->debug());
         this->_recvSection = new recv_section_type(this->comm(), this->_sendSection->getTag(), this->debug());
         this->_offsets     = new int[this->commSize()+1];
         this->_offsets[0]  = 0;
+        this->_localSize   = 0;
+      };
+    public:
+      Numbering(const Obj<topology_type>& topology, const std::string& label, int value) : ParallelObject(topology->comm(), topology->debug()), _topology(topology), _label(label), _value(value) {
+        this->__init();
+      };
+      template<typename Atlas>
+      Numbering(const Obj<Atlas>& atlas, const std::string& label, int value) : ParallelObject(atlas->comm(), atlas->debug()), _topology(atlas->getTopology()), _label(label), _value(value) {
+        this->__init();
       };
       ~Numbering() {
         delete [] this->_offsets;
     public:
       template<typename Atlas, typename Numbering>
       static Obj<Numbering> createIndices(const Obj<Atlas>& atlas, const Obj<Numbering>& numbering) {
-        Obj<Numbering> globalOffsets = new Numbering(numbering->getTopology(), numbering->getLabel(), numbering->getValue());
+        Obj<Numbering> globalOffsets = new Numbering(new Atlas(numbering->getTopology()), numbering->getLabel(), numbering->getValue());
         typename Atlas::patch_type patch = 0;
         // FIX: I think we can just generalize Numbering to take a stride based upon an Atlas
         //        However, then we also want a lightweight way of creating one numbering from another I think (maybe constructor?)
         const Obj<typename Numbering::topology_type::label_sequence>& points = numbering->getTopology()->getLabelStratum(patch, numbering->getLabel(), numbering->getValue());
         int offset = 0;
 
+        globalOffsets->copyCommunication(numbering);
+        globalOffsets->constructLocalOrder();
         for(typename Numbering::topology_type::label_sequence::iterator l_iter = points->begin(); l_iter != points->end(); ++l_iter) {
           if (numbering->isLocal(*l_iter)) {
             globalOffsets->setIndex(*l_iter, offset);
         globalOffsets->setLocalSize(offset);
         globalOffsets->calculateOffsets();
         globalOffsets->updateOrder();
-        globalOffsets->copyCommunication(numbering);
-        globalOffsets->fillSection();
-        globalOffsets->communicate();
-        globalOffsets->fillOrder();
+        //globalOffsets->fillSection();
+        //globalOffsets->communicate();
+        //globalOffsets->fillOrder();
+        globalOffsets->complete();
         return globalOffsets;
       };
     };

src/dm/mesh/sieve/Distribution.hh

       };
       #undef __FUNCT__
       #define __FUNCT__ "createScatterOverlap"
+      // This is the overlap for initially partitioned things where everyone communicates
+      //   We also need 0 --> everybody
+      //   An initial phase to figure out the selective overlap
       static void createScatterOverlap(const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
         const int rank = sendOverlap->commRank();
 
-        std::cout << "Rank: " << rank << " Size: " << sendOverlap->commSize() << std::endl;
         for(int p = 0; p < sendOverlap->commSize(); p++) {
           if (p != rank) {
-            // There are arrows to each rank whose color is the partition point (also the rank)
+            // Send arrow:   local point --- (remote point) ---> remote rank 
+            //   There are arrows to each rank whose color is the partition point (also the rank)
             sendOverlap->addCone(p, p, p);
-            //  The arrow is from rank 0 with partition point 0
-            recvOverlap->addCone(p, rank, rank);
+            // Receive arrow:   remote rank --- (remote point) ---> local point
+            recvOverlap->addCone(p, p, p);
           }
         }
         if (sendOverlap->debug) {sendOverlap->view(std::cout, "Send overlap for sieve scatter");}
             int c = 0;
 
             for(int p = 0; p < size; p++) {
-              sieve->addArrow(points[p], *b_iter, c++);
+              //sieve->addArrow(points[p], *b_iter, c++);
+              sieve->addArrow(points[p], *b_iter, c);
             }
           }
         }
-        sieve->stratify();
-        topology->stratify();
       }
       #undef __FUNCT__
-      #define __FUNCT__ "scatterSieve"
+      #define __FUNCT__ "updateOverlap"
+      template<typename SendSection, typename RecvSection>
+      static void updateOverlap(const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
+        const typename SendSection::topology_type::sheaf_type& sendRanks = sendSection->getAtlas()->getTopology()->getPatches();
+        const typename RecvSection::topology_type::sheaf_type& recvRanks = recvSection->getAtlas()->getTopology()->getPatches();
+
+        for(typename SendSection::topology_type::sheaf_type::const_iterator p_iter = sendRanks.begin(); p_iter != sendRanks.end(); ++p_iter) {
+          int                                                                       rank = p_iter->first;
+          const Obj<typename SendSection::topology_type::sieve_type::baseSequence>& base = p_iter->second->base();
+
+          for(typename SendSection::topology_type::sieve_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
+            const typename SendSection::value_type *points = sendSection->restrict(rank, *b_iter);
+            int size = sendSection->getAtlas()->getFiberDimension(rank, *b_iter);
+
+            for(int p = 0; p < size; p++) {
+              sendOverlap->addArrow(points[p], rank, points[p]);
+            }
+          }
+        }
+        for(typename RecvSection::topology_type::sheaf_type::const_iterator p_iter = recvRanks.begin(); p_iter != recvRanks.end(); ++p_iter) {
+          int                                                                       rank = p_iter->first;
+          const Obj<typename RecvSection::topology_type::sieve_type::baseSequence>& base = p_iter->second->base();
+
+          for(typename RecvSection::topology_type::sieve_type::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
+            const typename RecvSection::value_type *points = recvSection->restrict(rank, *b_iter);
+            int size = recvSection->getAtlas()->getFiberDimension(rank, *b_iter);
+
+            for(int p = 0; p < size; p++) {
+              recvOverlap->addArrow(rank, points[p], points[p]);
+            }
+          }
+        }
+      }
+      #undef __FUNCT__
+      #define __FUNCT__ "scatterTopology"
       static void scatterTopology(const Obj<topology_type>& oldTopology, const int numElements, const short assignment[], Obj<topology_type>& newTopology) {
         if (oldTopology->commSize() == 1) {
           newTopology = oldTopology;
         typedef typename ALE::New::OverlapValues<send_overlap_type, atlas_type, value_type> send_section_type;
         typedef typename ALE::New::OverlapValues<recv_overlap_type, atlas_type, value_type> recv_section_type;
         const Obj<send_section_type> sendSection = new send_section_type(oldTopology->comm(), oldTopology->debug());
-        const Obj<recv_section_type> recvSection = new recv_section_type(oldTopology->comm(), oldTopology->debug());
+        const Obj<recv_section_type> recvSection = new recv_section_type(oldTopology->comm(), sendSection->getTag(), oldTopology->debug());
 
         // Create Overlap - Could be calculated from partition, but would need communication
         createScatterOverlap(sendOverlap, recvOverlap);
         sectionCompletion::completeSection(sendOverlap, recvOverlap, coneSizeSection, coneSection, sendSection, recvSection);
         // Update sieve
         updateSieve(recvSection, newTopology);
+        newSieve->stratify();
+        newTopology->stratify();
+      };
+      static void coneCompletion(const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Obj<topology_type>& topology) {
+        typedef typename Mesh::point_type                                                   value_type;
+        typedef typename ALE::New::Completion<typename Mesh::topology_type, value_type>     completion_type;
+        typedef typename completion_type::atlas_type                                        atlas_type;
+        typedef typename ALE::New::OverlapValues<send_overlap_type, atlas_type, value_type> send_section_type;
+        typedef typename ALE::New::OverlapValues<recv_overlap_type, atlas_type, value_type> recv_section_type;
+        const Obj<send_section_type> sendSection = new send_section_type(topology->comm(), topology->debug());
+        const Obj<recv_section_type> recvSection = new recv_section_type(topology->comm(), sendSection->getTag(), topology->debug());
+
+        coneCompletion(sendOverlap, recvOverlap, topology, sendSection, recvSection);
+      };
+      #undef __FUNCT__
+      #define __FUNCT__ "coneCompletion"
+      template<typename SendSection, typename RecvSection>
+      static void coneCompletion(const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Obj<topology_type>& topology, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
+        if (sendOverlap->commSize() == 1) return;
+        typedef typename ALE::New::Completion<typename Mesh::topology_type, typename Mesh::point_type> completion_type;
+
+        // Distribute cones
+        const typename topology_type::patch_type               patch           = 0;
+        const Obj<completion_type::topology_type>              secTopology     = completion_type::createSendTopology(sendOverlap);
+        const Obj<typename completion_type::cone_size_section> coneSizeSection = new typename completion_type::cone_size_section(secTopology, topology->getPatch(patch));
+        const Obj<typename completion_type::cone_section>      coneSection     = new typename completion_type::cone_section(secTopology, topology->getPatch(patch));
+        sectionCompletion::completeSection(sendOverlap, recvOverlap, coneSizeSection, coneSection, sendSection, recvSection);
+        // Update cones
+        updateSieve(recvSection, topology);
       };
     };
   }

src/dm/mesh/sieve/Mesh.hh

     typedef ALE::New::Atlas<topology_type, ALE::Point>  atlas_type;
     typedef ALE::New::Section<atlas_type, double>       section_type;
     typedef std::map<std::string, Obj<section_type> >   SectionContainer;
-    typedef ALE::New::Numbering<topology_type>          numbering_type;
+    typedef ALE::New::NewNumbering<atlas_type>          numbering_type;
     typedef std::map<int, Obj<numbering_type> >         NumberingContainer;
     typedef std::map<std::string, Obj<numbering_type> > OrderContainer;
     typedef ALE::New::Section<atlas_type, ALE::pair<int,double> > foliated_section_type;
     typedef ALE::New::Section<atlas_type, ALE::pair<point_type, split_value> > split_section_type;
     typedef ALE::New::Completion<topology_type, point_type>::send_overlap_type send_overlap_type;
     typedef ALE::New::Completion<topology_type, point_type>::recv_overlap_type recv_overlap_type;
+    typedef ALE::New::Completion<topology_type, point_type>::atlas_type        comp_atlas_type;
+    typedef ALE::New::OverlapValues<send_overlap_type, comp_atlas_type, point_type> send_section_type;
+    typedef ALE::New::OverlapValues<recv_overlap_type, comp_atlas_type, point_type> recv_section_type;
     int debug;
   private:
     Obj<sieve_type>            topology;
     };
     const Obj<numbering_type>& getNumbering(const int depth) {
       if (this->numberings.find(depth) == this->numberings.end()) {
-        Obj<numbering_type> numbering = new numbering_type(this->getTopologyNew(), "depth", depth);
+        Obj<numbering_type> numbering = new numbering_type(new atlas_type(this->getTopologyNew()), "depth", depth);
         numbering->construct();
 
         std::cout << "Creating new numbering: " << depth << std::endl;
     };
     const Obj<numbering_type>& getLocalNumbering(const int depth) {
       if (this->localNumberings.find(depth) == this->localNumberings.end()) {
-        Obj<numbering_type> numbering = new numbering_type(this->getTopologyNew(), "depth", depth);
+        Obj<numbering_type> numbering = new numbering_type(new atlas_type(this->getTopologyNew()), "depth", depth);
         numbering->constructLocalOrder();
 
         std::cout << "Creating new local numbering: " << depth << std::endl;

src/dm/mesh/sieve/Sifter.hh

     template<class sourceInputSequence> 
     void 
     addCone(const Obj<sourceInputSequence>& sources, const typename traits::target_type& target, const typename traits::color_type& color){
-      if (debug) {std::cout << "Adding a cone " << std::endl;}
+      if (debug > 1) {std::cout << "Adding a cone " << std::endl;}
       for(typename sourceInputSequence::iterator iter = sources->begin(); iter != sources->end(); ++iter) {
-        if (debug) {std::cout << "Adding arrow from " << *iter << " to " << target << "(" << color << ")" << std::endl;}
+        if (debug > 1) {std::cout << "Adding arrow from " << *iter << " to " << target << "(" << color << ")" << std::endl;}
         this->addArrow(*iter, target, color);
       }
     };