dune-grid-glue  2.4.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 namespace Dune {
32 namespace 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 Dune::conditional<(side==0), typename Glue::Grid0View, typename Glue::Grid1View>::type GridView;
52  typedef typename Dune::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 Dune::conditional<(side==0), typename Glue::Grid0View, typename Glue::Grid1View>::type GridView;
186  typedef typename Dune::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 
275 public:
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
Write remote intersections to a vtk file for debugging purposes.
Definition: gridgluevtkwriter.hh:36
Definition: extractor.hh:54
static void write(const Glue &glue, const std::string &filenameTrunk)
Definition: gridgluevtkwriter.hh:277
Definition: gridglue.hh:33
Provides codimension-independent methods for grid extraction.
Definition: extractor.hh:47
Definition: extractor.hh:53