Dune Core Modules (2.4.1)

printgrid.hh
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
13
14namespace 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
Base class for Dune-Exceptions.
Definition: exceptions.hh:91
vector space out of a tensor product of fields.
Definition: fvector.hh:94
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:25
Default exception class for I/O errors.
Definition: exceptions.hh:256
Mesh entities of codimension 0 ("elements") allow to visit all intersections with "neighboring" eleme...
Definition: intersectioniterator.hh:84
Multiple codim and multiple geometry type mapper for leaf entities.
Definition: mcmgmapper.hh:297
A real mpi helper.
Definition: mpihelper.hh:160
int rank() const
return rank of process
Definition: mpihelper.hh:227
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:243
I trunc(const T &val, typename EpsilonType< T >::Type epsilon)
truncate using epsilon
Definition: float_cmp.cc:381
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:132
Mapper for multiple codim and multiple geometry types.
Helpers for dealing with MPI.
Dune namespace.
Definition: alignment.hh:10
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
Static tag representing a codimension.
Definition: dimension.hh:22
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 22, 23:33, 2024)