DUNE PDELab (git)

printgrid.hh
1// SPDX-FileCopyrightText: Copyright © 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
14
15namespace Dune {
16
17 namespace {
18
19 template<int dim>
20 struct ElementDataLayout
21 {
23 {
24 return gt.dim()==dim;
25 }
26 };
27
28 template<int dim>
29 struct NodeDataLayout
30 {
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
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
Default exception class for I/O errors.
Definition: exceptions.hh:231
Mesh entities of codimension 0 ("elements") allow to visit all intersections with "neighboring" eleme...
Definition: intersectioniterator.hh:83
A real mpi helper.
Definition: mpihelper.hh:181
int rank() const
return rank of process
Definition: mpihelper.hh:294
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
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
I trunc(const T &val, typename EpsilonType< T >::Type epsilon)
truncate using epsilon
Definition: float_cmp.cc:407
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:158
MCMGLayout mcmgElementLayout()
layout for elements (codim-0 entities)
Definition: mcmgmapper.hh:97
MCMGLayout mcmgVertexLayout()
layout for vertices (dim-0 entities)
Definition: mcmgmapper.hh:107
Mapper for multiple codim and multiple geometry types.
Helpers for dealing with MPI.
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
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
constexpr std::bool_constant<((II==value)||...)> contains(std::integer_sequence< T, II... >, std::integral_constant< T, value >)
Checks whether or not a given sequence contains a value.
Definition: integersequence.hh:137
Static tag representing a codimension.
Definition: dimension.hh:24
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)