dune-grid  2.4.1
printgrid.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 #ifndef DUNE_PRINTGRID_HH
4 #define DUNE_PRINTGRID_HH
5 
6 #include <fstream>
7 #include <string>
8 #include <sstream>
9 
10 #include <dune/common/exceptions.hh>
11 #include <dune/common/parallel/mpihelper.hh>
13 
14 namespace Dune {
15 
16  namespace {
17 
18  template<int dim>
19  struct ElementDataLayout
20  {
21  bool contains (Dune::GeometryType gt)
22  {
23  return gt.dim()==dim;
24  }
25  };
26 
27  template<int dim>
28  struct NodeDataLayout
29  {
30  bool contains (Dune::GeometryType gt)
31  {
32  return gt.dim()==0;
33  }
34  };
35 
36  // Move a point closer to basegeo's center by factor scale (used for drawing relative to the element)
37  template <typename B, typename C>
38  C centrify (const B& basegeo, const C& coords, const double scale) {
39  C ret = coords;
40  ret -= basegeo.center();
41  ret *= scale;
42  ret += basegeo.center();
43  return ret;
44  }
45 
46  // Add a line to the plotfile from p1 to p2
47  template <typename Coord>
48  void draw_line (std::ofstream &plotfile, const Coord &p1, const Coord &p2, std::string options) {
49  plotfile << "set object poly from ";
50  plotfile << p1[0] << "," << p1[1] << " to ";
51  plotfile << p2[0] << "," << p2[1] << " to ";
52  plotfile << p1[0] << "," << p1[1];
53  plotfile << " " << options << std::endl;
54  }
55 
56  }
57 
71  template <typename GridType>
72  void printGrid (const GridType& grid, const Dune::MPIHelper& helper, std::string output_file = "printgrid",
73  int size = 2000, bool execute_plot = true, bool png = true, bool local_corner_indices = true,
74  bool local_intersection_indices = true, bool outer_normals = true)
75  {
76 
77  // Create output file
78  std::stringstream of_name_stream(output_file);
79  of_name_stream << "_" << helper.rank();
80  output_file = of_name_stream.str();
81  std::string plot_file_name = output_file + ".gnuplot";
82  std::ofstream plotfile (plot_file_name, std::ios::out | std::ios::trunc);
83  if (!plotfile.is_open()) {
84  DUNE_THROW(Dune::IOError, "Could not create plot file " << output_file << "!");
85  return;
86  }
87 
88  // Basic plot settings
89  plotfile << "set size ratio -1" << std::endl;
90  if (png) {
91  plotfile << "set terminal png size " << size << "," << size << std::endl;
92  plotfile << "set output '" << output_file << ".png'" << std::endl;
93  } else {
94  plotfile << "set terminal svg size " << size << "," << size << " enhanced background rgb 'white'" << std::endl;
95  plotfile << "set output '" << output_file << ".svg'" << std::endl;
96  }
97 
98  // Get GridView
99  typedef typename GridType::LeafGridView GV;
100  const GV gv = grid.leafGridView();
101 
102  // Create mappers used to retrieve indices
105  ElementMapper elementmapper(grid);
106  NodeMapper nodemapper(grid);
107 
108  // Create iterators
109  typedef typename GV::template Codim<0 >::Iterator LeafIterator;
110  typedef typename GV::IntersectionIterator IntersectionIterator;
111 
112  LeafIterator it = gv.template begin<0>();
113 
114  // Will contain min/max coordinates. Needed for scaling of the plot
115  Dune::FieldVector<typename GridType::ctype,2> max_coord (it->geometry().center()), min_coord (max_coord);
116 
117  // Iterate over elements
118  for (; it != gv.template end<0>(); ++it) {
119 
120  const typename GV::template Codim<0>::Entity& entity = *it;
121  auto geo = entity.geometry();
122 
123  // Plot element index
124  int element_id = elementmapper.index(entity);
125  plotfile << "set label at " << geo.center()[0] << "," << geo.center()[1] << " '"
126  << element_id << "' center" << std::endl;
127 
128  for (int i = 0; i < geo.corners(); ++i) {
129  // Plot corner indices
130  const int globalNodeNumber1 = nodemapper.subIndex(entity, i, 2);
131  auto labelpos = centrify (geo, geo.corner(i), 0.7);
132  plotfile << "set label at " << labelpos[0] << "," << labelpos[1] << " '" << globalNodeNumber1;
133  if (local_corner_indices)
134  plotfile << "(" << i << ")";
135  plotfile << "' center" << std::endl;
136 
137  // Adapt min / max coordinates
138  for (int dim = 0; dim < 2; ++dim) {
139  if (geo.corner(i)[dim] < min_coord[dim])
140  min_coord[dim] = geo.corner(i)[dim];
141  else if (geo.corner(i)[dim] > max_coord[dim])
142  max_coord[dim] = geo.corner(i)[dim];
143  }
144  }
145 
146  // Iterate over intersections
147  for (IntersectionIterator is = gv.ibegin(entity); is != gv.iend(entity); ++is) {
148 
149  const typename GV::Intersection& intersection = *is;
150  auto igeo = intersection.geometry();
151 
152  // Draw intersection line
153  draw_line (plotfile, igeo.corner(0), igeo.corner(1), "fs empty border 1");
154 
155  // Plot local intersection index
156  if (local_intersection_indices) {
157  auto label_pos = centrify (geo, igeo.center(), 0.8);
158  plotfile << "set label at " << label_pos[0] << "," << label_pos[1]
159  << " '" << intersection.indexInInside() << "' center" << std::endl;
160  }
161 
162  // Plot outer normal
163  if (outer_normals) {
164  auto intersection_pos = igeo.center();
165  auto normal = intersection.centerUnitOuterNormal();
166  normal *= 0.15 * igeo.volume();
167  auto normal_end = intersection_pos + normal;
168  plotfile << "set arrow from " << intersection_pos[0] << "," << intersection_pos[1]
169  << " to " << normal_end[0] << "," << normal_end[1] << " lt rgb \"gray\"" << std::endl;
170  }
171 
172  // Get corners for inner intersection representation
173  auto inner_corner1 = centrify (geo, igeo.corner(0), 0.5);
174  auto inner_corner2 = centrify (geo, igeo.corner(1), 0.5);
175 
176  // Thick line in case of boundary()
177  if (intersection.boundary())
178  draw_line (plotfile, inner_corner1, inner_corner2, "fs empty border 3 lw 4");
179 
180  // Thin line with color according to neighbor()
181  if (intersection.neighbor())
182  draw_line (plotfile, inner_corner1, inner_corner2, "fs empty border 2");
183  else
184  draw_line (plotfile, inner_corner1, inner_corner2, "fs empty border 1");
185  }
186 
187  }
188 
189  // Finish plot, pass extend of the grid
190  Dune::FieldVector<typename GridType::ctype,2> extend (max_coord - min_coord);
191 
192  extend *= 0.2;
193  min_coord -= extend;
194  max_coord += extend;
195  plotfile << "plot [" << min_coord[0] << ":" << max_coord[0] << "] [" << min_coord[1]
196  << ":" << max_coord[1] << "] NaN notitle" << std::endl;
197  plotfile.close();
198 
199  if (execute_plot) {
200  std::string cmd = "gnuplot -p '" + plot_file_name + "'";
201  if (std::system (cmd.c_str()) != 0)
202  DUNE_THROW(Dune::Exception,"Error running GNUPlot: " << cmd);
203  }
204  }
205 
206 }
207 
208 #endif // #ifndef DUNE_PRINTGRID_HH
GeometryType
Type representing VTK&#39;s entity geometry types.
Definition: common.hh:178
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:72
Include standard header files.
Definition: agrid.hh:59
Mapper for multiple codim and multiple geometry types.
Mesh entities of codimension 0 ("elements") allow to visit all intersections with "neighboring" eleme...
Definition: common/grid.hh:360
Multiple codim and multiple geometry type mapper for leaf entities.
Definition: mcmgmapper.hh:295