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