33 typedef typename GridView::Grid
Grid;
36 static const int dimGrid = GridView::dimension;
39 typedef typename GridView::IndexSet IndexSet;
40 typedef typename GridView::template Codim< 0 >::Iterator ElementIterator;
41 typedef typename GridView::IntersectionIterator IntersectionIterator;
42 typedef typename GridView::template Codim< dimGrid >::EntityPointer VertexPointer;
44 typedef typename ElementIterator :: Entity Element ;
45 typedef typename Element :: EntityPointer ElementPointer;
46 typedef typename Element :: EntitySeed ElementSeed ;
48 typedef typename IndexSet::IndexType Index;
50 typedef Dune::ReferenceElement< typename Grid::ctype, dimGrid > RefElement;
51 typedef Dune::ReferenceElements< typename Grid::ctype, dimGrid > RefElements;
66 template <
class LoadBalanceHandle>
67 std::string
write (
const std::string &fileName,
const LoadBalanceHandle &ldb,
68 const int size,
const int rank )
const;
69 template <
class LoadBalanceHandle>
70 void write (
const std::string &fileName,
const LoadBalanceHandle &ldb,
71 const int size )
const;
82 const IndexSet& indexSet,
83 const Dune::GeometryType& elementType,
84 const std::vector< Index >& vertexIndex,
85 std::ostream &gridout )
const
88 if( element.type() != elementType )
92 const size_t vxSize = element.subEntities( Element::dimension );
93 Index vertices[ vxSize ];
94 for(
size_t i = 0; i < vxSize; ++i )
95 vertices[ i ] = vertexIndex[ indexSet.subIndex( element, i,
dimGrid ) ];
97 gridout << vertices[ 0 ];
98 for(
size_t i = 1; i < vxSize; ++i )
99 gridout <<
" " << vertices[ i ];
100 gridout << std::endl;
106inline std::string
DGFWriter< GV >::write (
const std::string &fileName,
const LoadBalanceHandle &ldb,
const int size,
const int rank )
const
108 std::stringstream newName;
109 newName << fileName <<
"." << size <<
"." << rank;
110 std::ofstream gridout( newName.str().c_str() );
112 gridout.setf( std::ios_base::scientific, std::ios_base::floatfield );
113 gridout.precision( 16 );
115 const IndexSet &indexSet = gridView_.indexSet();
118 gridout <<
"DGF" << std::endl;
120 const Index vxSize = indexSet.size( dimGrid );
121 std::vector< Index > vertexIndex( vxSize, vxSize );
123 gridout <<
"%" <<
" Elements = " << indexSet.size( 0 ) <<
" | Vertices = " << vxSize << std::endl;
126 gridout << std::endl <<
"VERTEX" << std::endl;
127 Index vertexCount = 0;
130 std::vector< ElementSeed > elementSeeds;
132 size_t countElements = 0 ;
134 const ElementIterator end = gridView_.template end< 0 >();
135 ElementIterator it = gridView_.template begin< 0 >();
136 for( ; it != end; ++it)
138 const Element& element = *it ;
139 if( ldb(element)==rank )
142 elementSeeds.push_back( element.seed() ) ;
146 const int numCorners = element.subEntities( dimGrid );
147 for(
int i=0; i<numCorners; ++i )
149 const Index vxIndex = indexSet.subIndex( element, i, dimGrid );
150 assert( vxIndex < vxSize );
151 if( vertexIndex[ vxIndex ] == vxSize )
153 vertexIndex[ vxIndex ] = vertexCount++;
154 gridout << element.geometry().corner( i ) << std::endl;
160 gridout <<
"#" << std::endl;
163 Dune::GeometryType simplex( Dune::GeometryType::simplex, dimGrid );
165 typedef typename std::vector<ElementSeed> :: const_iterator iterator ;
166 const iterator end = elementSeeds.end();
169 if( indexSet.size( simplex ) > 0 )
172 gridout << std::endl <<
"SIMPLEX" << std::endl;
176 for( iterator it = elementSeeds.begin(); it != end ; ++ it )
179 const ElementPointer ep = gridView_.grid().entityPointer( *it );
181 writeElement( *ep, indexSet, simplex, vertexIndex, gridout );
184 gridout <<
"#" << std::endl;
189 Dune::GeometryType cube( Dune::GeometryType::cube, dimGrid );
192 if( indexSet.size( cube ) > 0 )
195 gridout << std::endl <<
"CUBE" << std::endl;
196 for( iterator it = elementSeeds.begin(); it != end ; ++ it )
199 const ElementPointer ep = gridView_.grid().entityPointer( *it );
201 writeElement( *ep, indexSet, cube, vertexIndex, gridout );
205 gridout <<
"#" << std::endl;
210 gridout << std::endl <<
"BOUNDARYSEGMENTS" << std::endl;
211 for( iterator it = elementSeeds.begin(); it != end ; ++ it )
214 const ElementPointer ep = gridView_.grid().entityPointer( *it );
215 const Element& element = *ep;
216 if( !element.hasBoundaryIntersections() )
219 const RefElement &refElement = RefElements::general( element.type() );
221 const IntersectionIterator iend = gridView_.iend( element ) ;
222 for( IntersectionIterator iit = gridView_.ibegin( element ); iit != iend; ++iit )
224 if( !iit->boundary() )
227 const int boundaryId = iit->boundaryId();
228 if( boundaryId <= 0 )
230 std::cerr <<
"Warning: Ignoring nonpositive boundary id: "
231 << boundaryId <<
"." << std::endl;
235 const int faceNumber = iit->indexInInside();
236 const unsigned int faceSize = refElement.size( faceNumber, 1, dimGrid );
237 std::vector< Index > vertices( faceSize );
238 for(
unsigned int i = 0; i < faceSize; ++i )
240 const int j = refElement.subEntity( faceNumber, 1, i, dimGrid );
241 vertices[ i ] = vertexIndex[ indexSet.subIndex( element, j, dimGrid ) ];
243 gridout << boundaryId <<
" " << vertices[ 0 ];
244 for(
unsigned int i = 1; i < faceSize; ++i )
245 gridout <<
" " << vertices[ i ];
246 gridout << std::endl;
249 gridout <<
"#" << std::endl << std::endl;
251 gridout << std::endl <<
"GLOBALVERTEXINDEX" << std::endl;
252 for( iterator it = elementSeeds.begin(); it != end ; ++ it )
255 const ElementPointer ep = gridView_.grid().entityPointer( *it );
256 const Element& element = *ep;
257 const int numCorners = element.subEntities( dimGrid );
258 for(
int i=0; i<numCorners; ++i )
260 const Index vxIndex = indexSet.subIndex( element, i, dimGrid );
261 assert( vxIndex < vxSize );
262 if( vertexIndex[ vxIndex ] != vxSize )
264 vertexIndex[ vxIndex ] = vxSize;
265 gridout << vxIndex << std::endl;
269 gridout <<
"#" << std::endl << std::endl;
270 return newName.str();