dune-grid  2.2.1
dgfwriter.hh
Go to the documentation of this file.
1 #ifndef DUNE_DGFWRITER_HH
2 #define DUNE_DGFWRITER_HH
3 
9 #include <fstream>
10 #include <vector>
11 
12 #include <dune/grid/common/grid.hh>
13 #include <dune/geometry/referenceelements.hh>
14 
15 namespace Dune
16 {
17 
27  template< class GV >
28  class DGFWriter
29  {
30  typedef DGFWriter< GV > This;
31 
32  public:
34  typedef GV GridView;
36  typedef typename GridView::Grid Grid;
37 
39  static const int dimGrid = GridView::dimension;
40 
41  private:
42  typedef typename GridView::IndexSet IndexSet;
43  typedef typename GridView::template Codim< 0 >::Iterator ElementIterator;
44  typedef typename GridView::IntersectionIterator IntersectionIterator;
45  typedef typename GridView::template Codim< dimGrid >::EntityPointer VertexPointer;
46 
47  typedef typename IndexSet::IndexType Index;
48 
49  typedef GenericReferenceElement< typename Grid::ctype, dimGrid > RefElement;
50  typedef GenericReferenceElements< typename Grid::ctype, dimGrid > RefElements;
51 
52  public:
57  DGFWriter ( const GridView &gridView )
58  : gridView_( gridView )
59  {}
60 
65  void write ( std::ostream &gridout ) const;
66 
71  void write ( const std::string &fileName ) const;
72 
73  private:
74  GridView gridView_;
75  };
76 
77 
78 
79  template< class GV >
80  inline void DGFWriter< GV >::write ( std::ostream &gridout ) const
81  {
82  // set the stream to full double precision
83  gridout.setf( std::ios_base::scientific, std::ios_base::floatfield );
84  gridout.precision( 16 );
85 
86  const IndexSet &indexSet = gridView_.indexSet();
87 
88  // write DGF header
89  gridout << "DGF" << std::endl;
90 
91  const Index vxSize = indexSet.size( dimGrid );
92  std::vector< Index > vertexIndex( vxSize, vxSize );
93 
94  gridout << "%" << " Elements = " << indexSet.size( 0 ) << " | Vertices = " << vxSize << std::endl;
95 
96  // write all vertices into the "vertex" block
97  gridout << std::endl << "VERTEX" << std::endl;
98  Index vertexCount = 0;
99  typedef typename ElementIterator :: Entity Element ;
100  const ElementIterator end = gridView_.template end< 0 >();
101  for( ElementIterator it = gridView_.template begin< 0 >(); it != end; ++it )
102  {
103  const Element& element = *it ;
104  const int numCorners = element.template count< dimGrid > ();
105  for( int i=0; i<numCorners; ++i )
106  {
107  const Index vxIndex = indexSet.subIndex( element, i, dimGrid );
108  assert( vxIndex < vxSize );
109  if( vertexIndex[ vxIndex ] == vxSize )
110  {
111  vertexIndex[ vxIndex ] = vertexCount++;
112  gridout << element.geometry().corner( i ) << std::endl;
113  }
114  }
115  }
116  gridout << "#" << std::endl;
117  if( vertexCount != vxSize )
118  DUNE_THROW( GridError, "Index set reports wrong number of vertices." );
119 
120  if( dimGrid > 1 )
121  {
122  // write all simplices to the "simplex" block
123  gridout << std::endl << "SIMPLEX" << std::endl;
124  for( ElementIterator it = gridView_.template begin< 0 >(); it != end; ++it )
125  {
126  const Element& element = *it ;
127  if( !element.type().isSimplex() )
128  continue;
129 
130  std::vector< Index > vertices( dimGrid+1 );
131  for( size_t i = 0; i < vertices.size(); ++i )
132  vertices[ i ] = vertexIndex[ indexSet.subIndex( element, i, dimGrid ) ];
133 
134  gridout << vertices[ 0 ];
135  for( size_t i = 1; i < vertices.size(); ++i )
136  gridout << " " << vertices[ i ];
137  gridout << std::endl;
138  }
139  gridout << "#" << std::endl;
140  }
141 
142  // write all cubes to the "cube" block
143  gridout << std::endl << "CUBE" << std::endl;
144  for( ElementIterator it = gridView_.template begin< 0 >(); it != end; ++it )
145  {
146  const Element& element = *it ;
147  if( !element.type().isCube() )
148  continue;
149 
150  std::vector< Index > vertices( 1 << dimGrid );
151  for( size_t i = 0; i < vertices.size(); ++i )
152  vertices[ i ] = vertexIndex[ indexSet.subIndex( element, i, dimGrid ) ];
153 
154  gridout << vertices[ 0 ];
155  for( size_t i = 1; i < vertices.size(); ++i )
156  gridout << " " << vertices[ i ];
157  gridout << std::endl;
158  }
159  gridout << "#" << std::endl;
160 
161  // write all boundaries to the "boundarysegments" block
162  gridout << std::endl << "BOUNDARYSEGMENTS" << std::endl;
163  for( ElementIterator it = gridView_.template begin< 0 >(); it != end; ++it )
164  {
165  const Element& element = *it ;
166  if( !it->hasBoundaryIntersections() )
167  continue;
168 
169  const RefElement &refElement = RefElements::general( element.type() );
170 
171  const IntersectionIterator iend = gridView_.iend( element ) ;
172  for( IntersectionIterator iit = gridView_.ibegin( element ); iit != iend; ++iit )
173  {
174  if( !iit->boundary() )
175  continue;
176 
177  const int boundaryId = iit->boundaryId();
178  if( boundaryId <= 0 )
179  {
180  std::cerr << "Warning: Ignoring nonpositive boundary id: "
181  << boundaryId << "." << std::endl;
182  continue;
183  }
184 
185  const int faceNumber = iit->indexInInside();
186  const unsigned int faceSize = refElement.size( faceNumber, 1, dimGrid );
187  std::vector< Index > vertices( faceSize );
188  for( unsigned int i = 0; i < faceSize; ++i )
189  {
190  const int j = refElement.subEntity( faceNumber, 1, i, dimGrid );
191  vertices[ i ] = vertexIndex[ indexSet.subIndex( element, j, dimGrid ) ];
192  }
193  gridout << boundaryId << " " << vertices[ 0 ];
194  for( unsigned int i = 1; i < faceSize; ++i )
195  gridout << " " << vertices[ i ];
196  gridout << std::endl;
197  }
198  }
199  gridout << "#" << std::endl;
200 
201  // write the name into the "gridparameter" block
202  gridout << std::endl << "GRIDPARAMETER" << std::endl;
203  // gridout << "NAME " << gridView_.grid().name() << std::endl;
204  gridout << "#" << std::endl;
205 
206  gridout << std::endl << "#" << std::endl;
207  }
208 
209 
210  template< class GV >
211  inline void DGFWriter< GV >::write ( const std::string &fileName ) const
212  {
213  std::ofstream gridout( fileName.c_str() );
214  write( gridout );
215  }
216 
217 }
218 
219 #endif // #ifndef DUNE_DGFWRITER_HH