dune-grid  2.2.0
dgfwriter.hh
Go to the documentation of this file.
00001 #ifndef DUNE_DGFWRITER_HH
00002 #define DUNE_DGFWRITER_HH
00003 
00009 #include <fstream>
00010 #include <vector>
00011 
00012 #include <dune/grid/common/grid.hh>
00013 #include <dune/geometry/referenceelements.hh>
00014 
00015 namespace Dune
00016 {
00017 
00027   template< class GV >
00028   class DGFWriter
00029   {
00030     typedef DGFWriter< GV > This;
00031 
00032   public:
00034     typedef GV GridView;
00036     typedef typename GridView::Grid Grid;
00037 
00039     static const int dimGrid = GridView::dimension;
00040 
00041   private:
00042     typedef typename GridView::IndexSet IndexSet;
00043     typedef typename GridView::template Codim< 0 >::Iterator ElementIterator;
00044     typedef typename GridView::IntersectionIterator IntersectionIterator;
00045     typedef typename GridView::template Codim< dimGrid >::EntityPointer   VertexPointer;
00046 
00047     typedef typename IndexSet::IndexType Index;
00048 
00049     typedef GenericReferenceElement< typename Grid::ctype, dimGrid > RefElement;
00050     typedef GenericReferenceElements< typename Grid::ctype, dimGrid > RefElements;
00051 
00052   public:
00057     DGFWriter ( const GridView &gridView )
00058     : gridView_( gridView )
00059     {}
00060 
00065     void write ( std::ostream &gridout ) const;
00066 
00071     void write ( const std::string &fileName ) const;
00072 
00073   private:
00074     GridView gridView_;
00075   };
00076 
00077 
00078 
00079   template< class GV >
00080   inline void DGFWriter< GV >::write ( std::ostream &gridout ) const
00081   {
00082     // set the stream to full double precision
00083     gridout.setf( std::ios_base::scientific, std::ios_base::floatfield );
00084     gridout.precision( 16 );
00085     
00086     const IndexSet &indexSet = gridView_.indexSet();
00087    
00088     // write DGF header
00089     gridout << "DGF" << std::endl;
00090 
00091     const Index vxSize = indexSet.size( dimGrid );
00092     std::vector< Index > vertexIndex( vxSize, vxSize );
00093 
00094     gridout << "%" << " Elements = " << indexSet.size( 0 ) << "  |  Vertices = " << vxSize << std::endl;
00095 
00096     // write all vertices into the "vertex" block
00097     gridout << std::endl << "VERTEX" << std::endl;
00098     Index vertexCount = 0;
00099     typedef typename ElementIterator :: Entity  Element ;
00100     const ElementIterator end = gridView_.template end< 0 >();
00101     for( ElementIterator it = gridView_.template begin< 0 >(); it != end; ++it )
00102     {
00103       const Element& element = *it ;
00104       const int numCorners = element.template count< dimGrid > ();
00105       for( int i=0; i<numCorners; ++i ) 
00106       {
00107         const Index vxIndex = indexSet.subIndex( element, i, dimGrid );
00108         assert( vxIndex < vxSize );
00109         if( vertexIndex[ vxIndex ] == vxSize ) 
00110         {
00111           vertexIndex[ vxIndex ] = vertexCount++;
00112           gridout << element.geometry().corner( i ) << std::endl;
00113         }
00114       }
00115     }
00116     gridout << "#" << std::endl;
00117     if( vertexCount != vxSize )
00118       DUNE_THROW( GridError, "Index set reports wrong number of vertices." );
00119 
00120     // write all simplices to the "simplex" block
00121     gridout << std::endl << "SIMPLEX" << std::endl;
00122     for( ElementIterator it = gridView_.template begin< 0 >(); it != end; ++it )
00123     {
00124       const Element& element = *it ;
00125       if( ! element.type().isSimplex() )
00126         continue;
00127 
00128       std::vector< Index > vertices( dimGrid+1 );
00129       for( size_t i = 0; i < vertices.size(); ++i )
00130         vertices[ i ] = vertexIndex[ indexSet.subIndex( element, i, dimGrid ) ];
00131       
00132       gridout << vertices[ 0 ];
00133       for( size_t i = 1; i < vertices.size(); ++i )
00134         gridout << " " << vertices[ i ];
00135       gridout << std::endl;
00136     }
00137     gridout << "#" << std::endl;
00138 
00139     // write all cubes to the "cube" block
00140     gridout << std::endl << "CUBE" << std::endl;
00141     for( ElementIterator it = gridView_.template begin< 0 >(); it != end; ++it )
00142     {
00143       const Element& element = *it ;
00144       if( !element.type().isCube() )
00145         continue;
00146 
00147       std::vector< Index > vertices( 1 << dimGrid );
00148       for( size_t i = 0; i < vertices.size(); ++i )
00149         vertices[ i ] = vertexIndex[ indexSet.subIndex( element, i, dimGrid ) ];
00150       
00151       gridout << vertices[ 0 ];
00152       for( size_t i = 1; i < vertices.size(); ++i )
00153         gridout << " " << vertices[ i ];
00154       gridout << std::endl;
00155     }
00156     gridout << "#" << std::endl;
00157 
00158     // write all boundaries to the "boundarysegments" block
00159     gridout << std::endl << "BOUNDARYSEGMENTS" << std::endl;
00160     for( ElementIterator it = gridView_.template begin< 0 >(); it != end; ++it )
00161     {
00162       const Element& element = *it ;
00163       if( !it->hasBoundaryIntersections() )
00164         continue;
00165 
00166       const RefElement &refElement = RefElements::general( element.type() );
00167 
00168       const IntersectionIterator iend = gridView_.iend( element ) ;
00169       for( IntersectionIterator iit = gridView_.ibegin( element ); iit != iend; ++iit )
00170       {
00171         if( !iit->boundary() )
00172           continue;
00173 
00174         const int boundaryId = iit->boundaryId();
00175         if( boundaryId <= 0 )
00176         {
00177           std::cerr << "Warning: Ignoring nonpositive boundary id: "
00178                     << boundaryId << "." << std::endl;
00179           continue;
00180         }
00181 
00182         const int faceNumber = iit->indexInInside();
00183         const unsigned int faceSize = refElement.size( faceNumber, 1, dimGrid );
00184         std::vector< Index > vertices( faceSize );
00185         for( unsigned int i = 0; i < faceSize; ++i )
00186         {
00187           const int j = refElement.subEntity( faceNumber, 1, i, dimGrid );
00188           vertices[ i ] = vertexIndex[ indexSet.subIndex( element, j, dimGrid ) ];
00189         }
00190         gridout << boundaryId << "   " << vertices[ 0 ];
00191         for( unsigned int i = 1; i < faceSize; ++i )
00192           gridout << " " << vertices[ i ];
00193         gridout << std::endl;
00194       }
00195     }
00196     gridout << "#" << std::endl;
00197 
00198     // write the name into the "gridparameter" block
00199     gridout << std::endl << "GRIDPARAMETER" << std::endl;
00200     // gridout << "NAME " << gridView_.grid().name() << std::endl;
00201     gridout << "#" << std::endl;
00202     
00203     gridout << std::endl << "#" << std::endl;
00204   }
00205 
00206 
00207   template< class GV >
00208   inline void DGFWriter< GV >::write ( const std::string &fileName ) const
00209   {
00210     std::ofstream gridout( fileName.c_str() );
00211     write( gridout );
00212   }
00213 
00214 }
00215 
00216 #endif // #ifndef DUNE_DGFWRITER_HH