dgfwriter.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_DGFWRITER_HH
#define DUNE_DGFWRITER_HH

#include <fstream>
#include <vector>

#include <dune/grid/common/grid.hh>
#include <dune/geometry/referenceelements.hh>

namespace Dune
{

  template< class GV >
  class DGFWriter
  {
    typedef DGFWriter< GV > This;

  public:
    typedef GV GridView;
    typedef typename GridView::Grid Grid;

    static const int dimGrid = GridView::dimension;

  private:
    typedef typename GridView::IndexSet IndexSet;
    typedef typename GridView::template Codim< 0 >::Iterator ElementIterator;
    typedef typename GridView::IntersectionIterator IntersectionIterator;

    typedef typename ElementIterator :: Entity Element ;
    typedef typename Element :: EntitySeed ElementSeed ;

    typedef typename IndexSet::IndexType Index;

    typedef ReferenceElement< typename Grid::ctype, dimGrid > RefElement;
    typedef ReferenceElements< typename Grid::ctype, dimGrid > RefElements;

  public:
    DGFWriter ( const GridView &gridView )
      : gridView_( gridView )
    {}

    void write ( std::ostream &gridout,
                 const std::vector< Index >& newElemOrder,
                 const std::stringstream& addParams = std::stringstream() ) const;

    void write ( std::ostream &gridout ) const;

    void write ( std::ostream &gridout,
                 const std::stringstream& addParams ) const;

    void write ( const std::string &fileName ) const;

  protected:
    GridView gridView_;

  protected:
    //  helper methods

    // write all elements of type elementType
    void writeAllElements( const std::vector<ElementSeed>& elementSeeds,
                           const IndexSet& indexSet,
                           const GeometryType& elementType,
                           const std::vector< Index >& vertexIndex,
                           std::ostream &gridout ) const
    {
      if (!elementSeeds.empty()) {
        // perform grid traversal based on new element ordering
        for (const auto& seed : elementSeeds) {
          const Element element = gridView_.grid().entity(seed);
          writeElement(element, indexSet, elementType, vertexIndex, gridout);
        }
      }
      else {
        // perform default grid traversal
        for (const auto& element : elements(gridView_))
          writeElement(element, indexSet, elementType, vertexIndex, gridout);
      }
    }

    // write one element
    void writeElement( const Element& element,
                       const IndexSet& indexSet,
                       const GeometryType& elementType,
                       const std::vector< Index >& vertexIndex,
                       std::ostream &gridout ) const
    {
      // if element's type is not the same as the type to write the return
      if( element.type() != elementType )
        return ;

      // get vertex numbers of the element
      const size_t vxSize = element.subEntities( Element::dimension );
      std::vector<Index> vertices(vxSize);
      for( size_t i = 0; i < vxSize; ++i )
        vertices[ i ] = vertexIndex[ indexSet.subIndex( element, i, dimGrid ) ];

      gridout << vertices[ 0 ];
      for( size_t i = 1; i < vxSize; ++i )
        gridout << " " << vertices[ i ];
      gridout << std::endl;
    }
  };


  template< class GV >
  inline void DGFWriter< GV >::
  write ( std::ostream &gridout,
          const std::vector< Index >& newElemOrder,
          const std::stringstream& addParams ) const
  {
    // set the stream to full double precision
    gridout.setf( std::ios_base::scientific, std::ios_base::floatfield );
    gridout.precision( 16 );

    const IndexSet &indexSet = gridView_.indexSet();

    // vector containing entity seed (only needed if new ordering is given)
    std::vector< ElementSeed > elementSeeds;

    // if ordering was provided
    const size_t orderSize = newElemOrder.size() ;
    if( orderSize == indexSet.size( 0 ) )
    {
      const ElementIterator end = gridView_.template end< 0 >();
      ElementIterator it = gridView_.template begin< 0 >();

      if( it != end )
      {
        elementSeeds.resize( orderSize, (*it).seed() ) ;
        size_t countElements = 0 ;
        for( ; it != end; ++it, ++countElements )
        {
          const Element& element = *it ;
          assert( newElemOrder[ indexSet.index( element ) ] < orderSize );
          elementSeeds[ newElemOrder[ indexSet.index( element ) ] ] = element.seed();
        }

        // make sure that the size of the index set is equal
        // to the number of counted elements
        if( countElements != orderSize )
          DUNE_THROW(InvalidStateException,"DGFWriter::write: IndexSet not consecutive");
      }
    }

    // write DGF header
    gridout << "DGF" << std::endl;

    const Index vxSize = indexSet.size( dimGrid );
    std::vector< Index > vertexIndex( vxSize, vxSize );

    gridout << "%" << " Elements = " << indexSet.size( 0 ) << "  |  Vertices = " << vxSize << std::endl;

    // write all vertices into the "vertex" block
    gridout << std::endl << "VERTEX" << std::endl;
    Index vertexCount = 0;
    const ElementIterator end = gridView_.template end< 0 >();
    for( ElementIterator it = gridView_.template begin< 0 >(); it != end; ++it )
    {
      const Element& element = *it ;
      const int numCorners = element.subEntities( dimGrid );
      for( int i=0; i<numCorners; ++i )
      {
        const Index vxIndex = indexSet.subIndex( element, i, dimGrid );
        assert( vxIndex < vxSize );
        if( vertexIndex[ vxIndex ] == vxSize )
        {
          vertexIndex[ vxIndex ] = vertexCount++;
          gridout << element.geometry().corner( i ) << std::endl;
        }
      }
    }
    gridout << "#" << std::endl;
    if( vertexCount != vxSize )
      DUNE_THROW( GridError, "Index set reports wrong number of vertices." );

    if( dimGrid > 1 )
    {
      // type of element to write
      GeometryType simplex( GeometryType::simplex, dimGrid );

      // only write simplex block if grid view contains simplices
      if( indexSet.size( simplex ) > 0 )
      {
        // write all simplices to the "simplex" block
        gridout << std::endl << "SIMPLEX" << std::endl;

        // write all simplex elements
        writeAllElements( elementSeeds, indexSet, simplex, vertexIndex, gridout );

        // write end marker for block
        gridout << "#" << std::endl;
      }
    }

    {
      // cube geometry type
      GeometryType cube( GeometryType::cube, dimGrid );

      // only write cube block if grid view contains cubes
      if( indexSet.size( cube ) > 0 )
&#