dune-grid-glue 2.8-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
32
33namespace Dune {
34namespace 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
272public:
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