[bugfix][SFC] insert elements following the SFC ordering.

dune/alugrid/common/hsfc.hh

@@ -90,7 +90,7 @@ namespace ALUGridSFC {
 #endif // if HAVE_ZOLTAN
 
   template< class GridView >
-  void printSpaceFillingCurve ( const GridView& view, std::string name = "sfc" )
+  void printSpaceFillingCurve ( const GridView& view, std::string name = "sfc", const bool vtk = false  )
   {
     typedef typename GridView :: template Codim< 0 > :: Iterator  Iterator ;
     typedef typename Iterator :: Entity :: Geometry :: GlobalCoordinate GlobalCoordinate ;
@@ -120,35 +120,38 @@ namespace ALUGridSFC {
       gnuFile.close();
     }
 
-    std::stringstream vtkfilename;
-    vtkfilename << name << ".vtk";
-    if( view.grid().comm().size() > 1 )
-      vtkfilename << "." << view.grid().comm().rank();
-    std::ofstream vtkFile ( vtkfilename.str() );
-
-    if( vtkFile )
+    if( vtk )
     {
-      vtkFile << "# vtk DataFile Version 1.0" << std::endl;
-      vtkFile << "Line representation of vtk" << std::endl;
-      vtkFile << "ASCII" << std::endl;
-      vtkFile << "DATASET POLYDATA" << std::endl;
-      vtkFile << "POINTS "<< vertices.size() << " FLOAT" << std::endl;
+      std::stringstream vtkfilename;
+      vtkfilename << name << ".vtk";
+      if( view.grid().comm().size() > 1 )
+        vtkfilename << "." << view.grid().comm().rank();
+      std::ofstream vtkFile ( vtkfilename.str() );
 
-      for( size_t i=0; i<vertices.size(); ++i )
+      if( vtkFile )
       {
-        vtkFile << vertices[ i ];
-        for( int d=GlobalCoordinate::dimension; d<3; ++d )
-          vtkFile << " 0";
-        vtkFile << std::endl;
+        vtkFile << "# vtk DataFile Version 1.0" << std::endl;
+        vtkFile << "Line representation of vtk" << std::endl;
+        vtkFile << "ASCII" << std::endl;
+        vtkFile << "DATASET POLYDATA" << std::endl;
+        vtkFile << "POINTS "<< vertices.size() << " FLOAT" << std::endl;
+
+        for( size_t i=0; i<vertices.size(); ++i )
+        {
+          vtkFile << vertices[ i ];
+          for( int d=GlobalCoordinate::dimension; d<3; ++d )
+            vtkFile << " 0";
+          vtkFile << std::endl;
+        }
+
+        // lines, #lines, #entries
+        vtkFile << "LINES " << vertices.size()-1 << " " << (vertices.size()-1)*3 << std::endl;
+
+        for( size_t i=0; i<vertices.size()-1; ++i )
+          vtkFile << "2 " << i << " " << i+1 << std::endl;
+
+        vtkFile.close();
       }
-
-      // lines, #lines, #entries
-      vtkFile << "LINES " << vertices.size()-1 << " " << (vertices.size()-1)*3 << std::endl;
-
-      for( size_t i=0; i<vertices.size()-1; ++i )
-        vtkFile << "2 " << i << " " << i+1 << std::endl;
-
-      vtkFile.close();
     }
   }
 

dune/alugrid/common/structuredgridfactory.hh

@@ -102,6 +102,20 @@ namespace Dune
       }
 
     public:
+      template< class Entity >
+      double index( const Entity &entity ) const
+      {
+        alugrid_assert ( Entity::codimension == 0 );
+#ifdef USE_ALUGRID_SFC_ORDERING
+        // get center of entity's geometry
+        VertexType center = entity.geometry().center();
+        // get hilbert index in [0,1]
+        return sfc_.index( center );
+#else
+        return double(indexSet_.index( entity ));
+#endif
+      }
+
       template< class Entity >
       int rank( const Entity &entity ) const
       {
@@ -131,7 +145,6 @@ namespace Dune
         return pSize_-1;
       }
 
-    protected:
       void calculateElementCuts()
       {
         const size_t nElements = indexSet_.size( 0 );
@@ -196,7 +209,7 @@ namespace Dune
       std::vector< long int > elementCuts_ ;
 
 #ifdef USE_ALUGRID_SFC_ORDERING
-      // get element to hilbert index mapping
+      // get element to hilbert (or Z) index mapping
       SpaceFillingCurveOrdering< VertexType > sfc_;
 #endif
       const double maxIndex_ ;
@@ -371,19 +384,31 @@ namespace Dune
 
       // map global vertex ids to local ones
       std::map< IndexType, unsigned int > vtxMap;
+      std::map< double, const Entity > sortedElementList;
 
       const int numVertices = (1 << dim);
       std::vector< unsigned int > vertices( numVertices );
 
-      int nextElementIndex = 0;
       const ElementIterator end = gridView.template end< 0 >();
       for( ElementIterator it = gridView.template begin< 0 >(); it != end; ++it )
       {
         const Entity &entity = *it;
+
         // if the element does not belong to our partition, continue
         if( partitioner.rank( entity ) != myrank )
           continue;
 
+        const double elIndex = partitioner.index( entity );
+        assert( sortedElementList.find( elIndex ) == sortedElementList.end() );
+        sortedElementList.insert( std::make_pair( elIndex, entity ) );
+      }
+
+      int nextElementIndex = 0;
+      const auto endSorted = sortedElementList.end();
+      for( auto it = sortedElementList.begin(); it != endSorted; ++it )
+      {
+        const Entity &entity = (*it).second;
+
         // insert vertices and element
         const typename Entity::Geometry geo = entity.geometry();
         alugrid_assert( numVertices == geo.corners() );
@@ -397,6 +422,7 @@ namespace Dune
             factory.insertVertex( geo.corner( i ), vtxId );
           vertices[ i ] = result.first->second;
         }
+
         factory.insertElement( entity.type(), vertices );
         const int elementIndex = nextElementIndex++;