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