18 #ifndef DUNE_GRIDGLUE_ADAPTER_GRIDGLUEVTKWRITER_HH
19 #define DUNE_GRIDGLUE_ADAPTER_GRIDGLUEVTKWRITER_HH
24 #include <type_traits>
27 #include <dune/common/classname.hh>
28 #include <dune/geometry/type.hh>
29 #include <dune/geometry/referenceelements.hh>
44 template <
class Glue,
int s
ide>
45 static void writeExtractedPart(
const Glue& glue,
const std::string& filename)
47 static_assert((side==0 || side==1),
"'side' can only be 0 or 1");
51 fgrid.open(filename.c_str());
53 using GridView =
typename Glue::template GridView<side>;
54 using Extractor =
typename Glue::template GridPatch<side>;
56 typedef typename GridView::ctype ctype;
58 const int domdimw = GridView::dimensionworld;
62 std::string coordinatePadding;
63 for (
int i=domdimw; i<3; i++)
64 coordinatePadding +=
" 0";
66 fgrid <<
"# vtk DataFile Version 2.0\nFilename: " << filename <<
"\nASCII" << std::endl;
70 std::vector<typename Extractor::Coords> coords;
71 glue.template patch<side>().getCoords(coords);
73 fgrid << ((patchDim==3) ?
"DATASET UNSTRUCTURED_GRID" :
"DATASET POLYDATA") << std::endl;
74 fgrid <<
"POINTS " << coords.size() <<
" " << Dune::className<ctype>() << std::endl;
76 for (
size_t i=0; i<coords.size(); i++)
77 fgrid << coords[i] << coordinatePadding << std::endl;
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);
89 unsigned int faceCornerCount = 0;
90 for (
size_t i=0; i<faces.size(); i++)
91 faceCornerCount += faces[i].size();
93 fgrid << ((patchDim==3) ?
"CELLS " :
"POLYGONS ")
94 << geometryTypes.size() <<
" " << geometryTypes.size() + faceCornerCount << std::endl;
96 for (
size_t i=0; i<faces.size(); i++) {
98 fgrid << faces[i].size();
102 if (geometryTypes[i].isSimplex()) {
103 for (
int j=0; j<patchDim+1; j++)
104 fgrid <<
" " << faces[i][j];
106 }
else if (geometryTypes[i].isQuadrilateral()) {
107 fgrid <<
" " << faces[i][0] <<
" " << faces[i][1]
108 <<
" " << faces[i][3] <<
" " << faces[i][2];
110 }
else if (geometryTypes[i].isPyramid()) {
111 fgrid <<
" " << faces[i][0] <<
" " << faces[i][1]
112 <<
" " << faces[i][3] <<
" " << faces[i][2] <<
" " << faces[i][4];
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];
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];
125 DUNE_THROW(Dune::NotImplemented,
"Geometry type " << geometryTypes[i] <<
" not supported yet");
136 fgrid <<
"CELL_TYPES " << geometryTypes.size() << std::endl;
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;
148 DUNE_THROW(Dune::NotImplemented,
"Geometry type " << geometryTypes[i] <<
" not supported yet");
157 ctype accum = 0.0, delta = 1.0 / (ctype) (gridSubEntityData.size()-1);
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;
163 for (
typename GridSubEntityData::const_iterator sEIt = gridSubEntityData.begin();
164 sEIt != gridSubEntityData.end();
165 ++sEIt, accum += delta)
168 fgrid << accum << std::endl;
178 template <
class Glue,
int s
ide>
179 static void writeIntersections(
const Glue& glue,
const std::string& filename)
181 static_assert((side==0 || side==1),
"'side' can only be 0 or 1");
183 std::ofstream fmerged;
185 fmerged.open(filename.c_str());
187 using GridView =
typename Glue::template GridView<side>;
188 typedef typename GridView::ctype ctype;
190 const int domdimw = GridView::dimensionworld;
191 const int intersectionDim = Glue::Intersection::mydim;
194 std::string coordinatePadding;
195 for (
int i=domdimw; i<3; i++)
196 coordinatePadding +=
" 0";
198 int overlaps = glue.size();
202 using Extractor =
typename Glue::template GridPatch<0>;
203 std::vector<typename Extractor::Coords> coords;
204 glue.template patch<side>().
getCoords(coords);
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;
213 const auto& geometry = intersection.geometry();
214 for (
int i = 0; i < geometry.corners(); ++i)
215 fmerged << geometry.corner(i) << coordinatePadding << std::endl;
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);
226 unsigned int faceCornerCount = 0;
227 for (
size_t i=0; i<faces.size(); i++)
228 faceCornerCount += faces[i].size();
230 int grid0SimplexCorners = intersectionDim+1;
231 fmerged << ((intersectionDim==3) ?
"CELLS " :
"POLYGONS ")
232 << overlaps <<
" " << (grid0SimplexCorners+1)*overlaps << std::endl;
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;
242 if (intersectionDim==3) {
244 fmerged <<
"CELL_TYPES " << overlaps << std::endl;
246 for (
int i = 0; i < overlaps; i++)
247 fmerged <<
"10" << std::endl;
254 ctype accum = 0.0, delta = 1.0 / (ctype) (gridSubEntityData.size()-1);
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;
260 for (
typename GridSubEntityData::const_iterator sEIt = gridSubEntityData.begin();
261 sEIt != gridSubEntityData.end();
262 ++sEIt, accum += delta)
265 for (
int j = 0; j < sEIt->first.second; ++j)
266 fmerged << accum << std::endl;
273 template<
typename Glue>
274 static void write(
const Glue& glue,
const std::string& filenameTrunk)
278 writeExtractedPart<Glue,0>(glue,
279 filenameTrunk +
"-grid0.vtk");
281 writeIntersections<Glue,0>(glue,
282 filenameTrunk +
"-intersections-grid0.vtk");
285 writeExtractedPart<Glue,1>(glue,
286 filenameTrunk +
"-grid1.vtk");
288 writeIntersections<Glue,1>(glue,
289 filenameTrunk +
"-intersections-grid1.vtk");
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