dune-grid-glue 2.5-git
gridgluevtkwriter.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3/*
4 * Filename: GridGlueVtkWriter.hh
5 * Version: 1.0
6 * Created on: Mar 5, 2009
7 * Author: Gerrit Buse
8 * ---------------------------------
9 * Project: dune-grid-glue
10 * Description: Class thought to make graphical debugging of couplings easier.
11 *
12 */
18#ifndef DUNE_GRIDGLUE_ADAPTER_GRIDGLUEVTKWRITER_HH
19#define DUNE_GRIDGLUE_ADAPTER_GRIDGLUEVTKWRITER_HH
20
21
22#include <fstream>
23#include <iomanip>
24#include <type_traits>
25#include <vector>
26
27#include <dune/common/classname.hh>
28#include <dune/geometry/type.hh>
29#include <dune/geometry/referenceelements.hh>
30
31namespace Dune {
32namespace GridGlue {
33
37{
38
42 template <class Glue, int side>
43 static void writeExtractedPart(const Glue& glue, const std::string& filename)
44 {
45 static_assert((side==0 || side==1), "'side' can only be 0 or 1");
46
47 std::ofstream fgrid;
48
49 fgrid.open(filename.c_str());
50
51 typedef typename std::conditional<(side==0), typename Glue::Grid0View, typename Glue::Grid1View>::type GridView;
52 typedef typename std::conditional<(side==0), typename Glue::Grid0Patch, typename Glue::Grid1Patch>::type Extractor;
53
54 typedef typename GridView::ctype ctype;
55
56 const int domdimw = GridView::dimensionworld;
57 const int patchDim = Extractor::dim - Extractor::codim;
58
59 // coordinates have to be in R^3 in the VTK format
60 std::string coordinatePadding;
61 for (int i=domdimw; i<3; i++)
62 coordinatePadding += " 0";
63
64 fgrid << "# vtk DataFile Version 2.0\nFilename: " << filename << "\nASCII" << std::endl;
65
66 // WRITE POINTS
67 // ----------------
68 std::vector<typename Extractor::Coords> coords;
69 glue.template patch<side>().getCoords(coords);
70
71 fgrid << ((patchDim==3) ? "DATASET UNSTRUCTURED_GRID" : "DATASET POLYDATA") << std::endl;
72 fgrid << "POINTS " << coords.size() << " " << Dune::className<ctype>() << std::endl;
73
74 for (size_t i=0; i<coords.size(); i++)
75 fgrid << coords[i] << coordinatePadding << std::endl;
76
77 fgrid << std::endl;
78
79 // WRITE POLYGONS
80 // ----------------
81
82 std::vector<typename Extractor::VertexVector> faces;
83 std::vector<Dune::GeometryType> geometryTypes;
84 glue.template patch<side>().getFaces(faces);
85 glue.template patch<side>().getGeometryTypes(geometryTypes);
86
87 unsigned int faceCornerCount = 0;
88 for (size_t i=0; i<faces.size(); i++)
89 faceCornerCount += faces[i].size();
90
91 fgrid << ((patchDim==3) ? "CELLS " : "POLYGONS ")
92 << geometryTypes.size() << " " << geometryTypes.size() + faceCornerCount << std::endl;
93
94 for (size_t i=0; i<faces.size(); i++) {
95
96 fgrid << faces[i].size();
97
98 // vtk expects the vertices to by cyclically ordered
99 // therefore unfortunately we have to deal with several element types on a case-by-case basis
100 if (geometryTypes[i].isSimplex()) {
101 for (int j=0; j<patchDim+1; j++)
102 fgrid << " " << faces[i][j];
103
104 } else if (geometryTypes[i].isQuadrilateral()) {
105 fgrid << " " << faces[i][0] << " " << faces[i][1]
106 << " " << faces[i][3] << " " << faces[i][2];
107
108 } else if (geometryTypes[i].isPyramid()) {
109 fgrid << " " << faces[i][0] << " " << faces[i][1]
110 << " " << faces[i][3] << " " << faces[i][2] << " " << faces[i][4];
111
112 } else if (geometryTypes[i].isPrism()) {
113 fgrid << " " << faces[i][0] << " " << faces[i][2] << " " << faces[i][1]
114 << " " << faces[i][3] << " " << faces[i][5] << " " << faces[i][4];
115
116 } else if (geometryTypes[i].isHexahedron()) {
117 fgrid << " " << faces[i][0] << " " << faces[i][1]
118 << " " << faces[i][3] << " " << faces[i][2]
119 << " " << faces[i][4] << " " << faces[i][5]
120 << " " << faces[i][7] << " " << faces[i][6];
121
122 } else {
123 DUNE_THROW(Dune::NotImplemented, "Geometry type " << geometryTypes[i] << " not supported yet");
124 }
125
126 fgrid << std::endl;
127 }
128
129 fgrid << std::endl;
130
131 // 3d VTK files need an extra section specifying the CELL_TYPES aka GeometryTypes
132 if (patchDim==3) {
133
134 fgrid << "CELL_TYPES " << geometryTypes.size() << std::endl;
135
136 for (size_t i=0; i<geometryTypes.size(); i++) {
137 if (geometryTypes[i].isSimplex())
138 fgrid << "10" << std::endl;
139 else if (geometryTypes[i].isHexahedron())
140 fgrid << "12" << std::endl;
141 else if (geometryTypes[i].isPrism())
142 fgrid << "13" << std::endl;
143 else if (geometryTypes[i].isPyramid())
144 fgrid << "14" << std::endl;
145 else
146 DUNE_THROW(Dune::NotImplemented, "Geometry type " << geometryTypes[i] << " not supported yet");
147
148 }
149
150 }
151
152#if 0
153 // WRITE CELL DATA
154 // ---------------
155 ctype accum = 0.0, delta = 1.0 / (ctype) (gridSubEntityData.size()-1);
156
157 fgrid << "CELL_DATA " << gridSubEntityData.size() << std::endl;
158 fgrid << "SCALARS property_coding " << Dune::className<ctype>() << " 1" << std::endl;
159 fgrid << "LOOKUP_TABLE default" << std::endl;
160
161 for (typename GridSubEntityData::const_iterator sEIt = gridSubEntityData.begin();
162 sEIt != gridSubEntityData.end();
163 ++sEIt, accum += delta)
164 {
165 // "encode" the parent with one color...
166 fgrid << accum << std::endl;
167 }
168#endif
169 fgrid.close();
170 }
171
172
176 template <class Glue, int side>
177 static void writeIntersections(const Glue& glue, const std::string& filename)
178 {
179 static_assert((side==0 || side==1), "'side' can only be 0 or 1");
180
181 std::ofstream fmerged;
182
183 fmerged.open(filename.c_str());
184
185 typedef typename std::conditional<(side==0), typename Glue::Grid0View, typename Glue::Grid1View>::type GridView;
186 typedef typename std::conditional<(side==0),
187 typename Glue::Grid0IntersectionIterator,
188 typename Glue::Grid1IntersectionIterator>::type RemoteIntersectionIterator;
189
190 typedef typename GridView::ctype ctype;
191
192 const int domdimw = GridView::dimensionworld;
193 const int intersectionDim = Glue::Intersection::mydim;
194
195 // coordinates have to be in R^3 in the VTK format
196 std::string coordinatePadding;
197 for (int i=domdimw; i<3; i++)
198 coordinatePadding += " 0";
199
200 int overlaps = glue.size();
201
202 // WRITE POINTS
203 // ----------------
204 typedef typename Glue::Grid0Patch Extractor;
205 std::vector<typename Extractor::Coords> coords;
206 glue.template patch<side>().getCoords(coords);
207
208 // the merged grid (i.e. the set of remote intersections
209 fmerged << "# vtk DataFile Version 2.0\nFilename: " << filename << "\nASCII" << std::endl;
210 fmerged << ((intersectionDim==3) ? "DATASET UNSTRUCTURED_GRID" : "DATASET POLYDATA") << std::endl;
211 fmerged << "POINTS " << overlaps*(intersectionDim+1) << " " << Dune::className<ctype>() << std::endl;
212
213 for (RemoteIntersectionIterator isIt = glue.template ibegin<side>();
214 isIt != glue.template iend<side>();
215 ++isIt)
216 {
217 for (int i = 0; i < isIt->geometry().corners(); ++i)
218 fmerged << isIt->geometry().corner(i) << coordinatePadding << std::endl;
219 }
220
221 // WRITE POLYGONS
222 // ----------------
223
224 std::vector<typename Extractor::VertexVector> faces;
225 std::vector<Dune::GeometryType> geometryTypes;
226 glue.template patch<side>().getFaces(faces);
227 glue.template patch<side>().getGeometryTypes(geometryTypes);
228
229 unsigned int faceCornerCount = 0;
230 for (size_t i=0; i<faces.size(); i++)
231 faceCornerCount += faces[i].size();
232
233 int grid0SimplexCorners = intersectionDim+1;
234 fmerged << ((intersectionDim==3) ? "CELLS " : "POLYGONS ")
235 << overlaps << " " << (grid0SimplexCorners+1)*overlaps << std::endl;
236
237 for (int i = 0; i < overlaps; ++i) {
238 fmerged << grid0SimplexCorners;
239 for (int j=0; j<grid0SimplexCorners; j++)
240 fmerged << " " << grid0SimplexCorners*i+j;
241 fmerged << std::endl;
242 }
243
244 // 3d VTK files need an extra section specifying the CELL_TYPES aka GeometryTypes
245 if (intersectionDim==3) {
246
247 fmerged << "CELL_TYPES " << overlaps << std::endl;
248
249 for (int i = 0; i < overlaps; i++)
250 fmerged << "10" << std::endl;
251
252 }
253
254#if 0
255 // WRITE CELL DATA
256 // ---------------
257 ctype accum = 0.0, delta = 1.0 / (ctype) (gridSubEntityData.size()-1);
258
259 fmerged << "CELL_DATA " << overlaps << std::endl;
260 fmerged << "SCALARS property_coding " << Dune::className<ctype>() << " 1" << std::endl;
261 fmerged << "LOOKUP_TABLE default" << std::endl;
262
263 for (typename GridSubEntityData::const_iterator sEIt = gridSubEntityData.begin();
264 sEIt != gridSubEntityData.end();
265 ++sEIt, accum += delta)
266 {
267 // ...and mark all of its merged grid parts with the same color
268 for (int j = 0; j < sEIt->first.second; ++j)
269 fmerged << accum << std::endl;
270 }
271#endif
272 fmerged.close();
273 }
274
275public:
276 template<typename Glue>
277 static void write(const Glue& glue, const std::string& filenameTrunk)
278 {
279
280 // Write extracted grid and remote intersection on the grid0-side
281 writeExtractedPart<Glue,0>(glue,
282 filenameTrunk + "-grid0.vtk");
283
284 writeIntersections<Glue,0>(glue,
285 filenameTrunk + "-intersections-grid0.vtk");
286
287 // Write extracted grid and remote intersection on the grid1-side
288 writeExtractedPart<Glue,1>(glue,
289 filenameTrunk + "-grid1.vtk");
290
291 writeIntersections<Glue,1>(glue,
292 filenameTrunk + "-intersections-grid1.vtk");
293
294 }
295
296};
297
298} /* namespace GridGlue */
299} /* namespace Dune */
300
301#endif // DUNE_GRIDGLUE_ADAPTER_GRIDGLUEVTKWRITER_HH
Definition gridglue.hh:30
Write remote intersections to a vtk file for debugging purposes.
Definition gridgluevtkwriter.hh:37
static void write(const Glue &glue, const std::string &filenameTrunk)
Definition gridgluevtkwriter.hh:277
Provides codimension-independent methods for grid extraction.
Definition extractor.hh:43
@ dim
Definition extractor.hh:48
@ codim
Definition extractor.hh:49
void getCoords(std::vector< Dune::FieldVector< ctype, dimworld > > &coords) const
getter for the coordinates array
Definition extractor.hh:270