diff --git a/dune/grid/common/mcmgmapper.hh b/dune/grid/common/mcmgmapper.hh
index db2976e1c968809df8760b63776be6a3dce10cc2..97a57f8c0271f9dd6b25d4040e87cb03a0406e64 100644
--- a/dune/grid/common/mcmgmapper.hh
+++ b/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)];
diff --git a/dune/grid/io/file/vtk/common.hh b/dune/grid/io/file/vtk/common.hh
index a184eea8bc47012530516080e488c6debf47b43e..9898095e037b7fb6f37f43756f9017211e81562d 100644
--- a/dune/grid/io/file/vtk/common.hh
+++ b/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);
diff --git a/dune/grid/io/file/vtk/vtkwriter.hh b/dune/grid/io/file/vtk/vtkwriter.hh
index 75e950e2a4c50c9c770f964c74c0ae795aec49ed..5d75a8607fbc100750a9268f0f5bb64c662f518f 100644
--- a/dune/grid/io/file/vtk/vtkwriter.hh
+++ b/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;
};