Merge branch 'feature/vtkwriter-polyhedron' into 'master'

Polyhedron cells for VTKWriter.

This feature implements the support for polyhedrons in VTKWriter. This MR depends on !103 .

See merge request !105

dune/grid/common/mcmgmapper.hh

@@ -159,7 +159,11 @@ namespace Dune
      */
     Index subIndex (const typename GV::template Codim<0>::Entity& e, int i, unsigned int codim) const
     {
-      GeometryType gt=ReferenceElements<double,GV::dimension>::general(e.type()).type(i,codim);
+      const GeometryType eType = e.type();
+      GeometryType gt = eType.isNone() ?
+        GeometryType( GeometryType::none, GV::dimension - codim ) :
+        ReferenceElements<double,GV::dimension>::general(eType).type(i,codim) ;
+      //GeometryType gt=ReferenceElements<double,GV::dimension>::general(e.type()).type(i,codim);
       assert(layout.contains(gt));
       return is.subIndex(e, i, codim) + offset[GlobalGeometryTypeIndex::index(gt)];
     }
@@ -205,7 +209,11 @@ namespace Dune
      */
     bool contains (const typename GV::template Codim<0>::Entity& e, int i, int cc, Index& result) const
     {
-      GeometryType gt=ReferenceElements<double,GV::dimension>::general(e.type()).type(i,cc);
+      const GeometryType eType = e.type();
+      const GeometryType gt = eType.isNone() ?
+        GeometryType( GeometryType::none, GV::dimension - cc ) :
+        ReferenceElements<double,GV::dimension>::general(eType).type(i,cc) ;
+      //GeometryType gt=ReferenceElements<double,GV::dimension>::general(e.type()).type(i,cc);
       if (not layout.contains(gt))
         return false;
       result = is.subIndex(e, i, cc) + offset[GlobalGeometryTypeIndex::index(gt)];

dune/grid/io/file/vtk/common.hh

@@ -184,7 +184,8 @@ namespace Dune
       tetrahedron = 10,
       hexahedron = 12,
       prism = 13,
-      pyramid = 14 //, polyhedron = 42
+      pyramid = 14,
+      polyhedron = 42
     };
 
     //! mapping from GeometryType to VTKGeometryType
@@ -207,7 +208,7 @@ namespace Dune
       if (t.isNone() )
       {
         if( t.dim() == 2 ) return polygon;
-        // if( t.dim() == 3 ) return polyhedron;
+        if( t.dim() == 3 ) return polyhedron;
       }
 
       DUNE_THROW(IOError,"VTKWriter: unsupported GeometryType " << t);

dune/grid/io/file/vtk/vtkwriter.hh

@@ -14,6 +14,7 @@
 
 #include <vector>
 #include <list>
+#include <map>
 
 #include <dune/common/typetraits.hh>
 #include <dune/common/exceptions.hh>
@@ -576,7 +577,8 @@ namespace Dune
     explicit VTKWriter ( const GridView &gridView,
                          VTK::DataMode dm = VTK::conforming )
       : gridView_( gridView ),
-        datamode( dm )
+        datamode( dm ),
+        polyhedralCellsPresent_( checkForPolyhedralCells() )
     { }
 
     /**
@@ -1262,19 +1264,173 @@ namespace Dune
       // types
       if (n>1)
       {
-        std::shared_ptr<VTK::DataArrayWriter<unsigned char> > p3
-          (writer.makeArrayWriter<unsigned char>("types", 1, ncells));
-        if(!p3->writeIsNoop())
-          for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
+        {
+          std::shared_ptr<VTK::DataArrayWriter<unsigned char> > p3
+            (writer.makeArrayWriter<unsigned char>("types", 1, ncells));
+
+          if(!p3->writeIsNoop())
           {
-            int vtktype = VTK::geometryType(it->type());
-            p3->write(vtktype);
+            for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
+            {
+              int vtktype = VTK::geometryType(it->type());
+              p3->write(vtktype);
+            }
           }
+        }
+
+
+        // if polyhedron cells found also cell faces need to be written
+        if( polyhedralCellsPresent_ )
+        {
+          writeCellFaces( writer );
+        }
       }
 
       writer.endCells();
     }
 
+  protected:
+    bool checkForPolyhedralCells() const
+    {
+      // check if polyhedron cells are present
+      for( const auto& geomType : gridView_.indexSet().types( 0 ) )
+      {
+        if( VTK::geometryType( geomType ) == VTK::polyhedron )
+        {
+          return true;
+        }
+      }
+      return false;
+    }
+
+    //! write the connectivity array
+    virtual void writeCellFaces(VTK::VTUWriter& writer)
+    {
+      if( ! faceVertices_ )
+      {
+        faceVertices_.reset( new std::pair< std::vector<int>, std::vector<int> > () );
+        // fill face vertex structure
+        fillFaceVertices( cornerBegin(), cornerEnd(), gridView_.indexSet(),
+                          faceVertices_->first, faceVertices_->second );
+      }
+
+      std::vector< int >& faces = faceVertices_->first;
+      std::vector< int >& faceOffsets = faceVertices_->second;
+      assert( int(faceOffsets.size()) == ncells );
+
+      {
+        std::shared_ptr<VTK::DataArrayWriter< int > > p4
+          (writer.makeArrayWriter< int >("faces", 1, faces.size()));
+        if(!p4->writeIsNoop())
+        {
+          for( const auto& face : faces )
+            p4->write( face );
+        }
+      }
+
+      {
+        std::shared_ptr<VTK::DataArrayWriter< int > > p5
+          (writer.makeArrayWriter< int >("faceoffsets", 1, ncells ));
+        if(!p5->writeIsNoop())
+        {
+          for( const auto& offset : faceOffsets )
+            p5->write( offset );
+
+          // clear face vertex structure
+          faceVertices_.reset();
+        }
+      }
+    }
+
+    template <class CornerIterator, class IndexSet, class T>
+    inline void fillFaceVertices( CornerIterator it,
+                           const CornerIterator end,
+                           const IndexSet& indexSet,
+                           std::vector<T>& faces,
+                           std::vector<T>& faceOffsets )
+    {
+      if( n == 3 && it != end )
+      {
+        // clear output arrays
+        faces.clear();
+        faces.reserve( 15 * ncells );
+        faceOffsets.clear();
+        faceOffsets.reserve( ncells );
+
+        int offset = 0;
+
+        Cell element = *it;
+        int elIndex = indexSet.index( element );
+        std::vector< T > vertices;
+        vertices.reserve( 30 );
+        for( ; it != end; ++it )
+        {
+          const Cell& cell = *it ;
+          const int cellIndex = indexSet.index( cell ) ;
+          if( elIndex != cellIndex )
+          {
+            fillFacesForElement( element, indexSet, vertices, offset, faces, faceOffsets );
+
+            vertices.clear();
+            element = cell ;
+            elIndex = cellIndex ;
+          }
+          vertices.push_back( it.id() );
+        }
+
+        // fill faces for last element
+        fillFacesForElement( element, indexSet, vertices, offset, faces, faceOffsets );
+      }
+    }
+
+    template <class Entity, class IndexSet, class T>
+    static void fillFacesForElement( const Entity& element,
+                                     const IndexSet& indexSet,
+                                     const std::vector<T>& vertices,
+                                     T& offset,
+                                     std::vector<T>& faces,
+                                     std::vector<T>& faceOffsets )
+    {
+      const int dim = n;
+
+      std::map< T, T > vxMap;
+
+      // get number of local faces
+      const int nVertices = element.subEntities( dim );
+      for( int vx = 0; vx < nVertices; ++ vx )
+      {
+        const int vxIdx = indexSet.subIndex( element, vx, dim );
+        vxMap[ vxIdx ] = vertices[ vx ];
+      }
+
+      // get number of local faces
+      const int nFaces = element.subEntities( 1 );
+      // store number of faces for current element
+      faces.push_back( nFaces );
+      ++offset;
+      // extract each face as a set of vertex indices
+      for( int fce = 0; fce < nFaces; ++ fce )
+      {
+        // obtain face
+        const auto face = element.template subEntity< 1 > ( fce );
+
+        // get all vertex indices from current face
+        const int nVxFace = face.subEntities( dim );
+        faces.push_back( nVxFace );
+        ++offset ;
+        for( int i=0; i<nVxFace; ++i )
+        {
+          const T vxIndex = indexSet.subIndex( face, i, dim );
+          assert( vxMap.find( vxIndex ) != vxMap.end() );
+          faces.push_back( vxMap[ vxIndex ] );
+          ++offset ;
+        }
+      }
+
+      // store face offset for each element
+      faceOffsets.push_back( offset );
+    }
+
   protected:
     // the list of registered functions
     std::list<VTKLocalFunction> celldata;
@@ -1293,6 +1449,13 @@ namespace Dune
     // hold its number in the iteration order (VertexIterator)
     std::vector<int> number;
     VTK::DataMode datamode;
+
+    // true if polyhedral cells are present in the grid
+    const bool polyhedralCellsPresent_;
+
+    // pointer holding face vertex connectivity if needed
+    std::shared_ptr< std::pair< std::vector<int>, std::vector<int> > > faceVertices_;
+
   protected:
     VTK::OutputType outputtype;
   };