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