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