3 #ifndef DUNE_PRINTGRID_HH 4 #define DUNE_PRINTGRID_HH 9 #include <dune/common/exceptions.hh> 10 #include <dune/common/parallel/mpihelper.hh> 18 struct ElementDataLayout
36 template <
typename B,
typename C>
37 C centrify (
const B& basegeo,
const C& coords,
const double scale) {
39 ret -= basegeo.center();
41 ret += basegeo.center();
46 template <
typename Coord>
47 void draw_line (std::ofstream &plotfile,
const Coord &p1,
const Coord &p2, std::string options) {
48 plotfile <<
"set object poly from ";
49 plotfile << p1[0] <<
"," << p1[1] <<
" to ";
50 plotfile << p2[0] <<
"," << p2[1] <<
" to ";
51 plotfile << p1[0] <<
"," << p1[1];
52 plotfile <<
" " << options << std::endl;
70 template <
typename Gr
idType>
71 void printGrid (
const GridType& grid,
const Dune::MPIHelper& helper, std::string output_file =
"printgrid",
72 int size = 2000,
bool execute_plot =
true,
bool png =
true,
bool local_corner_indices =
true,
73 bool local_intersection_indices =
true,
bool outer_normals =
true)
77 output_file = output_file +
"_" + std::to_string(helper.rank());
78 std::string plot_file_name = output_file +
".gnuplot";
79 std::ofstream plotfile (plot_file_name, std::ios::out | std::ios::trunc);
80 if (!plotfile.is_open()) {
81 DUNE_THROW(Dune::IOError,
"Could not create plot file " << output_file <<
"!");
86 plotfile <<
"set size ratio -1" << std::endl;
88 plotfile <<
"set terminal png size " << size <<
"," << size << std::endl;
89 plotfile <<
"set output '" << output_file <<
".png'" << std::endl;
91 plotfile <<
"set terminal svg size " << size <<
"," << size <<
" enhanced background rgb 'white'" << std::endl;
92 plotfile <<
"set output '" << output_file <<
".svg'" << std::endl;
96 typedef typename GridType::LeafGridView GV;
97 const GV gv = grid.leafGridView();
105 typedef typename GV::template Codim<0 >::Iterator LeafIterator;
108 LeafIterator it = gv.template begin<0>();
111 Dune::FieldVector<typename GridType::ctype,2> max_coord (it->geometry().center()), min_coord (max_coord);
114 for (; it != gv.template end<0>(); ++it) {
116 const auto& entity = *it;
117 auto geo = entity.geometry();
120 int element_id = elementmapper.
index(entity);
121 plotfile <<
"set label at " << geo.center()[0] <<
"," << geo.center()[1] <<
" '" 122 << element_id <<
"' center" << std::endl;
124 for (
int i = 0; i < geo.corners(); ++i) {
126 const int globalNodeNumber1 = nodemapper.
subIndex(entity, i, 2);
127 auto labelpos = centrify (geo, geo.corner(i), 0.7);
128 plotfile <<
"set label at " << labelpos[0] <<
"," << labelpos[1] <<
" '" << globalNodeNumber1;
129 if (local_corner_indices)
130 plotfile <<
"(" << i <<
")";
131 plotfile <<
"' center" << std::endl;
134 for (
int dim = 0; dim < 2; ++dim) {
135 if (geo.corner(i)[dim] < min_coord[dim])
136 min_coord[dim] = geo.corner(i)[dim];
137 else if (geo.corner(i)[dim] > max_coord[dim])
138 max_coord[dim] = geo.corner(i)[dim];
143 for (IntersectionIterator is = gv.ibegin(entity); is != gv.iend(entity); ++is) {
145 const auto& intersection = *is;
146 auto igeo = intersection.geometry();
149 draw_line (plotfile, igeo.corner(0), igeo.corner(1),
"fs empty border 1");
152 if (local_intersection_indices) {
153 auto label_pos = centrify (geo, igeo.center(), 0.8);
154 plotfile <<
"set label at " << label_pos[0] <<
"," << label_pos[1]
155 <<
" '" << intersection.indexInInside() <<
"' center" << std::endl;
160 auto intersection_pos = igeo.center();
161 auto normal = intersection.centerUnitOuterNormal();
162 normal *= 0.15 * igeo.volume();
163 auto normal_end = intersection_pos + normal;
164 plotfile <<
"set arrow from " << intersection_pos[0] <<
"," << intersection_pos[1]
165 <<
" to " << normal_end[0] <<
"," << normal_end[1] <<
" lt rgb \"gray\"" << std::endl;
169 auto inner_corner1 = centrify (geo, igeo.corner(0), 0.5);
170 auto inner_corner2 = centrify (geo, igeo.corner(1), 0.5);
173 if (intersection.boundary())
174 draw_line (plotfile, inner_corner1, inner_corner2,
"fs empty border 3 lw 4");
177 if (intersection.neighbor())
178 draw_line (plotfile, inner_corner1, inner_corner2,
"fs empty border 2");
180 draw_line (plotfile, inner_corner1, inner_corner2,
"fs empty border 1");
186 Dune::FieldVector<typename GridType::ctype,2> extend (max_coord - min_coord);
191 plotfile <<
"plot [" << min_coord[0] <<
":" << max_coord[0] <<
"] [" << min_coord[1]
192 <<
":" << max_coord[1] <<
"] NaN notitle" << std::endl;
196 std::string cmd =
"gnuplot -p '" + plot_file_name +
"'";
197 if (std::system (cmd.c_str()) != 0)
198 DUNE_THROW(Dune::Exception,
"Error running GNUPlot: " << cmd);
204 #endif // #ifndef DUNE_PRINTGRID_HH Index index(const EntityType &e) const
Map entity to array index.
Definition: mapper.hh:119
void printGrid(const GridType &grid, const Dune::MPIHelper &helper, std::string output_file="printgrid", int size=2000, bool execute_plot=true, bool png=true, bool local_corner_indices=true, bool local_intersection_indices=true, bool outer_normals=true)
Print a grid as a gnuplot for testing and development.
Definition: printgrid.hh:71
Index subIndex(const typename G::Traits::template Codim< 0 >::Entity &e, int i, unsigned int codim) const
Map subentity i of codim cc of a codim 0 entity to array index.
Definition: mapper.hh:133
MCMGLayout mcmgElementLayout()
layout for elements (codim-0 entities)
Definition: mcmgmapper.hh:150
Mesh entities of codimension 0 ("elements") allow to visit all intersections with "neighboring" eleme...
Definition: common/grid.hh:345
MCMGLayout mcmgVertexLayout()
layout for vertices (dim-0 entities)
Definition: mcmgmapper.hh:160
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:178
Mapper for multiple codim and multiple geometry types.
Mapper interface.
Definition: mapper.hh:107
Multiple codim and multiple geometry type mapper for leaf entities.
Definition: mcmgmapper.hh:471
Include standard header files.
Definition: agrid.hh:58