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