18 #ifndef DUNE_GRIDGLUE_ADAPTER_GRIDGLUEVTKWRITER_HH
19 #define DUNE_GRIDGLUE_ADAPTER_GRIDGLUEVTKWRITER_HH
26 #include <dune/common/classname.hh>
27 #include <dune/common/typetraits.hh>
28 #include <dune/geometry/type.hh>
29 #include <dune/geometry/referenceelements.hh>
40 template <
class Glue,
int s
ide>
41 static void writeExtractedPart(
const Glue& glue,
const std::string& filename)
43 static_assert((side==0 || side==1),
"'side' can only be 0 or 1");
47 fgrid.open(filename.c_str());
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;
57 typedef typename GridView::ctype ctype;
59 const int dim = GridView::dimension;
60 const int domdimw = GridView::dimensionworld;
61 const int patchDim = Extractor::dim - Extractor::codim;
64 std::string coordinatePadding;
65 for (
int i=domdimw; i<3; i++)
66 coordinatePadding +=
" 0";
68 fgrid <<
"# vtk DataFile Version 2.0\nFilename: " << filename <<
"\nASCII" << std::endl;
72 std::vector<typename Extractor::Coords> coords;
73 glue.template patch<side>().getCoords(coords);
75 fgrid << ((patchDim==3) ?
"DATASET UNSTRUCTURED_GRID" :
"DATASET POLYDATA") << std::endl;
76 fgrid <<
"POINTS " << coords.size() <<
" " << Dune::className<ctype>() << std::endl;
78 for (
size_t i=0; i<coords.size(); i++)
79 fgrid << coords[i] << coordinatePadding << std::endl;
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);
91 unsigned int faceCornerCount = 0;
92 for (
size_t i=0; i<faces.size(); i++)
93 faceCornerCount += faces[i].size();
95 fgrid << ((patchDim==3) ?
"CELLS " :
"POLYGONS ")
96 << geometryTypes.size() <<
" " << geometryTypes.size() + faceCornerCount << std::endl;
98 for (
size_t i=0; i<faces.size(); i++) {
100 fgrid << faces[i].size();
104 if (geometryTypes[i].isSimplex()) {
105 for (
int j=0; j<dim; j++)
106 fgrid <<
" " << faces[i][j];
108 }
else if (geometryTypes[i].isQuadrilateral()) {
109 fgrid <<
" " << faces[i][0] <<
" " << faces[i][1]
110 <<
" " << faces[i][3] <<
" " << faces[i][2];
112 }
else if (geometryTypes[i].isPyramid()) {
113 fgrid <<
" " << faces[i][0] <<
" " << faces[i][1]
114 <<
" " << faces[i][3] <<
" " << faces[i][2] <<
" " << faces[i][4];
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];
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];
127 DUNE_THROW(Dune::NotImplemented,
"Geometry type " << geometryTypes[i] <<
" not supported yet");
138 fgrid <<
"CELL_TYPES " << geometryTypes.size() << std::endl;
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;
150 DUNE_THROW(Dune::NotImplemented,
"Geometry type " << geometryTypes[i] <<
" not supported yet");
159 ctype accum = 0.0, delta = 1.0 / (ctype) (gridSubEntityData.size()-1);
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;
165 for (
typename GridSubEntityData::const_iterator sEIt = gridSubEntityData.begin();
166 sEIt != gridSubEntityData.end();
167 ++sEIt, accum += delta)
170 fgrid << accum << std::endl;
180 template <
class Glue,
int s
ide>
181 static void writeIntersections(
const Glue& glue,
const std::string& filename)
183 static_assert((side==0 || side==1),
"'side' can only be 0 or 1");
185 std::ofstream fmerged;
187 fmerged.open(filename.c_str());
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;
196 typedef typename GridView::ctype ctype;
198 const int dim = GridView::dimension;
199 const int domdimw = GridView::dimensionworld;
200 const int intersectionDim = Glue::Intersection::mydim;
203 std::string coordinatePadding;
204 for (
int i=domdimw; i<3; i++)
205 coordinatePadding +=
" 0";
207 int overlaps = glue.size();
211 typedef typename Glue::Grid0Patch Extractor;
212 std::vector<typename Extractor::Coords> coords;
213 glue.template patch<side>().getCoords(coords);
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;
220 for (RemoteIntersectionIterator isIt = glue.template ibegin<side>();
221 isIt != glue.template iend<side>();
224 for (
int i = 0; i < isIt->geometry().corners(); ++i)
225 fmerged << isIt->geometry().corner(i) << coordinatePadding << std::endl;
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);
236 unsigned int faceCornerCount = 0;
237 for (
size_t i=0; i<faces.size(); i++)
238 faceCornerCount += faces[i].size();
240 int grid0SimplexCorners = dim-Glue::Grid0Patch::codim+1;
241 fmerged << ((intersectionDim==3) ?
"CELLS " :
"POLYGONS ")
242 << overlaps <<
" " << (grid0SimplexCorners+1)*overlaps << std::endl;
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;
252 if (intersectionDim==3) {
254 fmerged <<
"CELL_TYPES " << overlaps << std::endl;
256 for (
size_t i=0; i<overlaps; i++)
257 fmerged <<
"10" << std::endl;
264 ctype accum = 0.0, delta = 1.0 / (ctype) (gridSubEntityData.size()-1);
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;
270 for (
typename GridSubEntityData::const_iterator sEIt = gridSubEntityData.begin();
271 sEIt != gridSubEntityData.end();
272 ++sEIt, accum += delta)
275 for (
int j = 0; j < sEIt->first.second; ++j)
276 fmerged << accum << std::endl;
283 template<
typename Glue>
284 static void write(
const Glue& glue,
const std::string& filenameTrunk)
288 writeExtractedPart<Glue,0>(glue,
289 filenameTrunk +
"-grid0.vtk");
291 writeIntersections<Glue,0>(glue,
292 filenameTrunk +
"-intersections-grid0.vtk");
295 writeExtractedPart<Glue,1>(glue,
296 filenameTrunk +
"-grid1.vtk");
298 writeIntersections<Glue,1>(glue,
299 filenameTrunk +
"-intersections-grid1.vtk");
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