Commits

Matt Knepley committed 8a198a9

Added the interpolation switch to overlapTest2
CoSieve.hh:
- removed atlasNew from UniformSection
- fixed restrictPoint() in UniformSection
- Added updateAdd() to UniformSection
Completion.hh:
- Added a new Numbering class

Hg-commit: 660d09fa9c3a8f02e05140ecc578e541f840df62

Comments (0)

Files changed (4)

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

 typedef ALE::New::NewNumbering<atlas_type>        new_numbering_type;
 
 typedef struct {
-  int debug; // The debugging level
+  int        debug;       // The debugging level
+  PetscTruth interpolate; // Construct missing elements of the mesh
 } Options;
 
 #undef __FUNCT__
 // This test has
 //  - A mesh overlapping itself
 //  - Single points overlapping single points
-PetscErrorCode DoubletTest(const Obj<send_section_type>& sendSection, const Obj<recv_section_type>& recvSection, Options *options)
+PetscErrorCode DoubletTest(MPI_Comm comm, Options *options)
 {
+  Obj<send_section_type> sendSection = new send_section_type(comm, options->debug);
+  Obj<recv_section_type> recvSection = new recv_section_type(comm, sendSection->getTag(), options->debug);
   Obj<topology_type>     topology    = new topology_type(sendSection->comm(), options->debug);
   Obj<send_overlap_type> sendOverlap = new send_overlap_type(sendSection->comm(), options->debug);
   Obj<recv_overlap_type> recvOverlap = new recv_overlap_type(recvSection->comm(), options->debug);
   Obj<numbering_type> numbering = new numbering_type(topology, "depth", 0);
 
   PetscFunctionBegin;
-  ALE::Test::OverlapTest::constructDoublet2(topology);
   numbering->construct();
   numbering->view("Vertex numbering");
   PetscFunctionReturn(0);
   Obj<new_numbering_type> numbering = new new_numbering_type(atlas, "depth", 0);
 
   PetscFunctionBegin;
-  ALE::Test::OverlapTest::constructDoublet2(topology);
   numbering->construct();
   numbering->view("Vertex numbering");
   PetscFunctionReturn(0);
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
-  options->debug = 0;
+  options->debug       = 0;
+  options->interpolate = PETSC_TRUE;
 
   ierr = PetscOptionsBegin(comm, "", "Options for sifter stress test", "Sieve");CHKERRQ(ierr);
     ierr = PetscOptionsInt("-debug", "The debugging level", "overlap1.c", options->debug, &options->debug, PETSC_NULL);CHKERRQ(ierr);
+    ierr = PetscOptionsTruth("-interpolate", "Construct missing elements of the mesh", "overlap1.c", options->interpolate, &options->interpolate, PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();
   PetscFunctionReturn(0);
 }
   comm = PETSC_COMM_WORLD;
   ierr = ProcessOptions(comm, &options);CHKERRQ(ierr);
   try {
-    Obj<send_section_type> sendSection = new send_section_type(comm, options.debug);
-    Obj<recv_section_type> recvSection = new recv_section_type(comm, sendSection->getTag(), options.debug);
-    Obj<topology_type>     topology    = new topology_type(comm, options.debug);
+    Obj<topology_type> topology = new topology_type(comm, options.debug);
 
-    //ierr = DoubletTest(sendSection, recvSection, &options);CHKERRQ(ierr);
-    ierr = NumberingTest(topology, &options);CHKERRQ(ierr);
+    ALE::Test::OverlapTest::constructDoublet2(topology, options.interpolate);
+    //ierr = DoubletTest(comm, &options);CHKERRQ(ierr);
     ierr = NewNumberingTest(topology, &options);CHKERRQ(ierr);
+    ierr = NumberingTest(topology, &options);CHKERRQ(ierr);
   } catch (ALE::Exception e) {
     std::cout << e << std::endl;
   }

src/dm/mesh/examples/tests/overlapTest.hh

       //       \ | /
       //        \|/
       //       1 | 1
-      static void constructDoublet2(const Obj<topology_type>& topology) {
+      static void constructDoublet2(const Obj<topology_type>& topology, bool interpolate = true) {
         Obj<sieve_type>     sieve = new sieve_type(topology->comm(), topology->debug());
         Obj<std::set<int> > cone  = new std::set<int>();
 
         if (topology->commRank() == 0) {
-          cone->insert(3);cone->insert(4);cone->insert(5);
-          sieve->addCone(cone, 8);cone->clear();
-          cone->insert(0);cone->insert(1);
-          sieve->addCone(cone, 3);cone->clear();
-          cone->insert(1);cone->insert(2);
-          sieve->addCone(cone, 4);cone->clear();
-          cone->insert(2);cone->insert(0);
-          sieve->addCone(cone, 5);cone->clear();
+          if (interpolate) {
+            cone->insert(3);cone->insert(4);cone->insert(5);
+            sieve->addCone(cone, 8);cone->clear();
+            cone->insert(0);cone->insert(1);
+            sieve->addCone(cone, 3);cone->clear();
+            cone->insert(1);cone->insert(2);
+            sieve->addCone(cone, 4);cone->clear();
+            cone->insert(2);cone->insert(0);
+            sieve->addCone(cone, 5);cone->clear();
+          } else {
+            cone->insert(0);cone->insert(1);cone->insert(2);
+            sieve->addCone(cone, 8);cone->clear();
+          }
         } else {
-          cone->insert(4);cone->insert(6);cone->insert(7);
-          sieve->addCone(cone, 9);cone->clear();
-          cone->insert(1);cone->insert(3);
-          sieve->addCone(cone, 6);cone->clear();
-          cone->insert(3);cone->insert(2);
-          sieve->addCone(cone, 7);cone->clear();
-          cone->insert(2);cone->insert(1);
-          sieve->addCone(cone, 4);cone->clear();
+          if (interpolate) {
+            cone->insert(4);cone->insert(6);cone->insert(7);
+            sieve->addCone(cone, 9);cone->clear();
+            cone->insert(1);cone->insert(3);
+            sieve->addCone(cone, 6);cone->clear();
+            cone->insert(3);cone->insert(2);
+            sieve->addCone(cone, 7);cone->clear();
+            cone->insert(2);cone->insert(1);
+            sieve->addCone(cone, 4);cone->clear();
+          } else {
+            cone->insert(1);cone->insert(2);cone->insert(3);
+            sieve->addCone(cone, 7);cone->clear();
+          }
         }
+        sieve->stratify();
         topology->setPatch(0, sieve);
         topology->stratify();
       };

src/dm/mesh/sieve/CoSieve.hh

       typedef std::map<patch_type, array_type>       values_type;
     protected:
       Obj<atlas_type> _atlas;
-      Obj<atlas_type> _atlasNew;
       values_type     _arrays;
     public:
       UniformSection(MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug) {
-        this->_atlas    = new atlas_type(comm, debug);
-        this->_atlasNew = NULL;
+        this->_atlas = new atlas_type(comm, debug);
       };
-      UniformSection(const Obj<topology_type>& topology) : ParallelObject(topology->comm(), topology->debug()), _atlasNew(NULL) {
+      UniformSection(const Obj<topology_type>& topology) : ParallelObject(topology->comm(), topology->debug()) {
         this->_atlas = new atlas_type(topology);
       };
-      UniformSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas), _atlasNew(NULL) {};
+      UniformSection(const Obj<atlas_type>& atlas) : ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas) {};
     protected:
       value_type *getRawArray(const int size) {
         static value_type *array   = NULL;
       };
       const value_type *restrictPoint(const patch_type& patch, const point_type& p) {
         this->checkPatch(patch);
-        return &this->_arrays[patch][p];
+        return this->_arrays[patch][p].v;
       };
       void update(const patch_type& patch, const point_type& p, const value_type v[]) {
         this->_atlas->checkPatch(patch);
 #endif
         }
       };
+      void updateAdd(const patch_type& patch, const point_type& p, const value_type v[]) {
+        this->_atlas->checkPatch(patch);
+        array_type& array = this->_arrays[patch];
+
+        if (this->getTopology()->height(patch) == 1) {
+          // Only avoids the copy of closure()
+          const int& dim = this->_atlas->restrict(patch, p)[0];
+          int        j   = -1;
+
+          for(int i = 0; i < dim; ++i) {
+            array[p].v[i] += v[++j];
+          }
+          // Should be closure()
+          const Obj<typename sieve_type::coneSequence>& cone = this->getTopology()->getPatch(patch)->cone(p);
+
+          for(typename sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
+            const int& dim = this->_atlas->restrict(patch, *p_iter)[0];
+
+            for(int i = 0; i < dim; ++i) {
+              array[*p_iter].v[i] += v[++j];
+            }
+          }
+        } else {
+          throw ALE::Exception("Not yet implemented for interpolated sieves");
+        }
+      };
       void updatePoint(const patch_type& patch, const point_type& p, const value_type v[]) {
         this->_atlas->checkPatch(patch);
         for(int i = 0; i < fiberDim; ++i) {

src/dm/mesh/sieve/Completion.hh

       };
     };
 
+    template<typename Topology_, typename Value_ = int>
+    class Numbering : public UniformSection<Topology_, Value_> {
+    public:
+      typedef UniformSection<Topology_, Value_>  base_type;
+      typedef Topology_                          topology_type;
+      typedef typename topology_type::patch_type patch_type;
+      typedef typename topology_type::sieve_type sieve_type;
+      typedef typename topology_type::point_type point_type;
+      typedef typename base_type::atlas_type     atlas_type;
+      typedef typename atlas_type::chart_type    chart_type;
+      typedef Value_                             value_type;
+      typedef typename base_type::array_type     array_type;
+      typedef typename base_type::values_type    values_type;
+      typedef typename ALE::Sifter<int,point_type,point_type> send_overlap_type;
+      typedef typename ALE::Sifter<point_type,int,point_type> recv_overlap_type;
+    protected:
+      std::string               _label;
+      value_type                _value;
+      Obj<send_overlap_type>    _sendOverlap;
+      Obj<recv_overlap_type>    _recvOverlap;
+      int                       _localSize;
+      int                      *_offsets;
+      std::map<int, point_type> _invOrder;
+    public:
+      Numbering(const Obj<topology_type>& topology, const std::string& label, int value) : UniformSection<Topology_, Value_>(topology), _label(label), _value(value), _localSize(0) {
+        this->_sendOverlap = new send_overlap_type(this->comm(), this->debug());
+        this->_recvOverlap = new recv_overlap_type(this->comm(), this->debug());
+        this->_offsets     = new int[this->commSize()+1];
+        this->_offsets[0]  = 0;
+      };
+      ~Numbering() {
+        delete [] this->_offsets;
+      };
+    public: // Accessors
+      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() const {return this->_sendOverlap;};
+      void setSendOverlap(const Obj<send_overlap_type>& overlap) {this->_sendOverlap = overlap;};
+      const Obj<recv_overlap_type>& getRecvOverlap() const {return this->_recvOverlap;};
+      void setRecvOverlap(const Obj<recv_overlap_type>& overlap) {this->_recvOverlap = overlap;};
+    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);};
+      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;};
+      point_type getPoint(const int& index) {return this->_invOrder[index];};
+      bool hasPoint(const point_type& point) {
+        const typename atlas_type::chart_type& chart = this->getChart(0);
+        return (chart->find(point) != chart->end());
+      };
+    public: // construction
+      void constructOverlap(const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
+        const patch_type patch = 0;
+        const Obj<typename topology_type::label_sequence>& points = this->getTopology()->getLabelStratum(patch, this->_label, this->_value);
+
+        point_type *sendBuf = new point_type[points->size()];
+        int         size    = 0;
+        for(typename topology_type::label_sequence::iterator l_iter = points->begin(); l_iter != points->end(); ++l_iter) {
+          sendBuf[size++] = *l_iter;
+        }
+        int *sizes   = new int[this->commSize()];
+        int *offsets = new int[this->commSize()+1];
+        point_type *remotePoints = NULL;
+        int        *remoteRanks  = NULL;
+
+        // Change to Allgather() for the correct binning algorithm
+        MPI_Gather(&size, 1, MPI_INT, sizes, 1, MPI_INT, 0, this->comm());
+        if (this->commRank() == 0) {
+          offsets[0] = 0;
+          for(int p = 1; p <= this->commSize(); p++) {
+            offsets[p] = offsets[p-1] + sizes[p-1];
+          }
+          remotePoints = new point_type[offsets[this->commSize()]];
+        }
+        MPI_Gatherv(sendBuf, size, MPI_INT, remotePoints, sizes, offsets, MPI_INT, 0, this->comm());
+        std::map<int, std::map<int, std::set<point_type> > > overlapInfo;
+
+        if (this->commRank() == 0) {
+          for(int p = 0; p < this->commSize(); p++) {
+            std::sort(&remotePoints[offsets[p]], &remotePoints[offsets[p+1]]);
+          }
+          for(int p = 0; p < this->commSize(); p++) {
+            for(int q = p+1; q < this->commSize(); q++) {
+              std::set_intersection(&remotePoints[offsets[p]], &remotePoints[offsets[p+1]],
+                                    &remotePoints[offsets[q]], &remotePoints[offsets[q+1]],
+                                    std::insert_iterator<std::set<point_type> >(overlapInfo[p][q], overlapInfo[p][q].begin()));
+              overlapInfo[q][p] = overlapInfo[p][q];
+            }
+            sizes[p]     = overlapInfo[p].size()*2;
+            offsets[p+1] = offsets[p] + sizes[p];
+          }
+          remoteRanks = new int[offsets[this->commSize()]];
+          int       k = 0;
+          for(int p = 0; p < this->commSize(); p++) {
+            for(typename std::map<int, std::set<point_type> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
+              remoteRanks[k*2]   = r_iter->first;
+              remoteRanks[k*2+1] = r_iter->second.size();
+              k++;
+            }
+          }
+        }
+        int numOverlaps;
+        MPI_Scatter(sizes, 1, MPI_INT, &numOverlaps, 1, MPI_INT, 0, this->comm());
+        int *overlapRanks = new int[numOverlaps];
+        MPI_Scatterv(remoteRanks, sizes, offsets, MPI_INT, overlapRanks, numOverlaps, MPI_INT, 0, this->comm());
+        if (this->commRank() == 0) {
+          for(int p = 0, k = 0; p < this->commSize(); p++) {
+            sizes[p] = 0;
+            for(int r = 0; r < (int) overlapInfo[p].size(); r++) {
+              sizes[p] += remoteRanks[k*2+1];
+              k++;
+            }
+            offsets[p+1] = offsets[p] + sizes[p];
+          }
+          for(int p = 0, k = 0; p < this->commSize(); p++) {
+            for(typename std::map<int, std::set<point_type> >::iterator r_iter = overlapInfo[p].begin(); r_iter != overlapInfo[p].end(); ++r_iter) {
+              int rank = r_iter->first;
+              for(typename std::set<point_type>::iterator p_iter = (overlapInfo[p][rank]).begin(); p_iter != (overlapInfo[p][rank]).end(); ++p_iter) {
+                remotePoints[k++] = *p_iter;
+              }
+            }
+          }
+        }
+        int numOverlapPoints = 0;
+        for(int r = 0; r < numOverlaps/2; r++) {
+          numOverlapPoints += overlapRanks[r*2+1];
+        }
+        point_type *overlapPoints = new point_type[numOverlapPoints];
+        MPI_Scatterv(remotePoints, sizes, offsets, MPI_INT, overlapPoints, numOverlapPoints, MPI_INT, 0, this->comm());
+
+        for(int r = 0, k = 0; r < numOverlaps/2; r++) {
+          int rank = overlapRanks[r*2];
+
+          for(int p = 0; p < overlapRanks[r*2+1]; p++) {
+            point_type point = overlapPoints[k++];
+
+            sendOverlap->addArrow(point, rank, point);
+            recvOverlap->addArrow(rank, point, point);
+          }
+        }
+
+        delete [] overlapPoints;
+        delete [] overlapRanks;
+        delete [] sizes;
+        delete [] offsets;
+        if (this->commRank() == 0) {
+          delete [] remoteRanks;
+          delete [] remotePoints;
+        }
+        if (this->debug()) {
+          sendOverlap->view("Send overlap");
+          recvOverlap->view("Receive overlap");
+        }
+      };
+      // Number all locals points with the given label value
+      //   points in the overlap are only numbered by the owner with the lowest rank
+      void constructLocalOrder(const Obj<send_overlap_type>& sendOverlap) {
+        const patch_type patch = 0;
+        const Obj<typename topology_type::label_sequence>& points = this->getTopology()->getLabelStratum(patch, this->_label, this->_value);
+
+        this->setFiberDimensionByLabel(patch, this->_label, this->_value, 1);
+        this->_localSize = 0;
+        for(typename topology_type::label_sequence::iterator l_iter = points->begin(); l_iter != points->end(); ++l_iter) {
+          int val;
+
+          if (sendOverlap->capContains(*l_iter)) {
+            const Obj<typename send_overlap_type::traits::supportSequence>& sendPatches = sendOverlap->support(*l_iter);
+            int minRank = sendOverlap->commSize();
+
+            for(typename send_overlap_type::traits::supportSequence::iterator p_iter = sendPatches->begin(); p_iter != sendPatches->end(); ++p_iter) {
+              if (*p_iter < minRank) minRank = *p_iter;
+            }
+            if (minRank < sendOverlap->commRank()) {
+              val = -1;
+            } else {
+              val = this->_localSize++;
+            }
+          } else {
+            val = this->_localSize++;
+          }
+          this->update(patch, *l_iter, &val);
+        }
+      };
+      // Construct the inverse map from numbers to points
+      //   If we really need this, then we should consider using a label
+      void constructInverseOrder() {
+        const chart_type& patch = this->getAtlas()->getPatch(0);
+
+        for(typename chart_type::iterator p_iter = patch.begin(); p_iter != patch.end(); ++p_iter) {
+          this->_invOrder[this->getIndex(*p_iter)] = *p_iter;
+        }
+      };
+      // Calculate process offsets
+      void calculateOffsets() {
+        MPI_Allgather(&this->_localSize, 1, MPI_INT, &(this->_offsets[1]), 1, MPI_INT, this->comm());
+        for(int p = 2; p <= this->commSize(); p++) {
+          this->_offsets[p] += this->_offsets[p-1];
+        }
+      };
+      // Update local offsets based upon process offsets
+      void updateOrder() {
+        const patch_type patch = 0;
+        const Obj<typename topology_type::label_sequence>& points = this->getAtlas()->getTopology()->getLabelStratum(patch, this->_label, this->_value);
+        const int val = this->_offsets[this->commRank()];
+
+        for(typename topology_type::label_sequence::iterator l_iter = points->begin(); l_iter != points->end(); ++l_iter) {
+          if (this->isLocal(*l_iter)) {
+            this->updateAdd(patch, *l_iter, &val);
+          }
+        }
+      };
+      // Communicate numbers in the overlap
+      void complete(const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, 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;
+        typedef typename ALE::New::ConstantSection<topology_type, int> constant_sizer;
+        const Obj<send_section_type> sendSection = new send_section_type(this->comm(), this->debug());
+        const Obj<recv_section_type> recvSection = new recv_section_type(this->comm(), sendSection->getTag(), this->debug());
+        const Obj<constant_sizer>    sizer       = new constant_sizer(this->comm(), 1, this->debug());
+        Obj<Numbering<topology_type> > filler(this);
+        (*filler.refCnt)++;
+
+        Completion<topology_type, int>::completeSection(sendOverlap, recvOverlap, sizer, filler, sendSection, recvSection);
+        Obj<typename recv_overlap_type::traits::baseSequence> recvPoints = recvOverlap->base();
+
+        for(typename recv_overlap_type::traits::baseSequence::iterator r_iter = recvPoints->begin(); r_iter != recvPoints->end(); ++r_iter) {
+          this->setFiberDimension(0, *r_iter, 1);
+        }
+        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 = 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) && !allowDuplicates) {
+                ostringstream msg;
+                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);
+              this->update(0, *r_iter, &val);
+            }
+          }
+        }
+      };
+      void construct() {
+        this->constructOverlap(this->_sendOverlap, this->_recvOverlap);
+        this->construct(this->_sendOverlap, this->_recvOverlap);
+      };
+      void construct(const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap) {
+        this->constructLocalOrder(sendOverlap);
+        this->calculateOffsets();
+        this->updateOrder();
+        this->complete(sendOverlap, recvOverlap);
+      };
+      void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
+        ostringstream txt;
+        int rank;
+
+        if (comm == MPI_COMM_NULL) {
+          comm = this->comm();
+          rank = this->commRank();
+        } else {
+          MPI_Comm_rank(comm, &rank);
+        }
+        if (name == "") {
+          if(rank == 0) {
+            txt << "viewing a Numbering" << std::endl;
+          }
+        } else {
+          if(rank == 0) {
+            txt << "viewing Numbering '" << name << "'" << std::endl;
+          }
+        }
+        for(typename values_type::const_iterator a_iter = this->_arrays.begin(); a_iter != this->_arrays.end(); ++a_iter) {
+          const patch_type  patch = a_iter->first;
+          array_type&       array = this->_arrays[patch];
+
+          txt << "[" << this->commRank() << "]: Patch " << patch << std::endl;
+          const typename atlas_type::chart_type& chart = this->_atlas->getPatch(patch);
+
+          for(typename atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
+            const typename atlas_type::point_type& p   = *p_iter;
+            const typename atlas_type::value_type  dim = this->_atlas->restrict(patch, p)[0];
+
+            if (dim != 0) {
+              txt << "[" << this->commRank() << "]:   " << p << " --> ";
+              for(int i = 0; i < dim; i++) {
+                if (array[p].v[i] >= 0) {
+                  txt << array[p].v[i];
+                } else {
+                  txt << -(array[p].v[i]+1) << " (global)";
+                }
+              }
+              txt << std::endl;
+            }
+          }
+        }
+        PetscSynchronizedPrintf(comm, txt.str().c_str());
+        PetscSynchronizedFlush(comm);
+      };
+    };
+
     template<typename Atlas_>
     class NewNumbering : public Section<Atlas_, int> {
     public:
       };
     };
 
+#if 0
     template<typename Topology_>
     class Numbering : public ParallelObject {
     public:
         PetscSynchronizedFlush(this->comm());
       };
     };
+#endif
 
     class GlobalOrder {
     public: