OpenLB 1.7
Loading...
Searching...
No Matches
gnuplotHeatMapWriter.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2017 Marc Haussmann
4 * E-mail contact: info@openlb.net
5 * The most recent release of OpenLB can be downloaded at
6 * <http://www.openlb.net/>
7 *
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License
10 * as published by the Free Software Foundation; either version 2
11 * of the License, or (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public
19 * License along with this program; if not, write to the Free
20 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
21 * Boston, MA 02110-1301, USA.
22*/
23
24#ifndef GNUPLOT_HEATMAP_WRITER_HH
25#define GNUPLOT_HEATMAP_WRITER_HH
26
27#include <fstream>
28#include <iostream>
29#include <unistd.h>
30
31#include "core/singleton.h"
32#include "io/fileName.hh"
35
36namespace olb {
37
38namespace heatmap {
39
40template <typename T>
41void write(BlockReduction3D2D<T>& blockReduction, int iT, const plotParam<T> param, const std::vector<T>& valueArea)
42{
43 detail::genericHeatMapInterface(blockReduction, blockReduction, iT, valueArea, param);
44}
45
46template <typename T>
47void write(BlockReduction2D2D<T>& blockReduction, int iT, const plotParam<T> param, const std::vector<T>& valueArea)
48{
49 detail::genericHeatMapInterface(blockReduction.getPlaneDiscretizationIn3D(), blockReduction, iT,
50 valueArea, param);
51}
52
53namespace detail {
54
55template <typename T>
56void genericHeatMapInterface(const HyperplaneLattice3D<T>& hyperPlane, BlockF2D<T>& blockData, int iT,
57 const std::vector<T>& valueArea, const plotParam<T>& plot)
58{
59
60
61 if ( blockData.getTargetDim() != 1 ) {
62 std::cout << "HeatMap Error: Functor targetDim is not 1. " << std::endl;
63 exit(-1);
64 }
65 else {
66 if ( singleton::mpi().getRank() == 0 ) {
68 param.plot = &plot;
69 param.blockData = &blockData;
70 param.hyperPlane = &hyperPlane;
72 param.quantityname = param.blockData->getName();
73 param.quantityname.erase (param.quantityname.begin(), param.quantityname.begin()+15);
74 param.quantityname.pop_back();
75 std::string name;
76 if (plot.name == "") {
77 name = param.quantityname;
78 }
79 else {
80 name = plot.name;
81 }
82 //std::stringstream fNameStream;
83 //fNameStream << param.dir << "data/"<< name << "iT" << std::setw(7) << std::setfill('0') << iT;
84
85 param.matrixPath = createFileName( param.dir + "data/", name, iT) + ".matrix";
86 param.csvPath = createFileName( param.dir + "data/", name, iT) + ".csv";
87 param.jpegPath = createFileName( param.dir, name, iT) + ".jpeg";
88 param.pngPath = createFileName( param.dir, name, iT) + ".png";
89 param.plotFilePath = createFileName( param.dir + "data/", name, iT) +".p";
90
91 for (std::size_t pos = 0; pos < param.plotFilePath.length(); pos++) {
92 if (param.plotFilePath.at(pos) == '(' || param.plotFilePath.at(pos) == ')') {
93 param.plotFilePath.replace(pos, 1, "");
94 }
95 }
96
97 param.nx = param.hyperPlane->getNx();
98 param.ny = param.hyperPlane->getNy();
99 param.spacing = param.hyperPlane->getPhysSpacing();
100 param.origin = param.hyperPlane->getPhysOrigin();
101 param.normal = crossProduct3D(param.hyperPlane->getVectorU(),param.hyperPlane->getVectorV());
102 param.zoomMin = plot.zoomOrigin;
103 param.zoomMax = plot.zoomOrigin + plot.zoomExtend;
104 param.iT = iT;
105 //aspect ratio declared here
106 param.aspect = param.nx/(double)param.ny;
107 if (plot.fullScreenPlot) {
108 T base = std::min(param.nx, param.ny);
109 param.canvasX = param.aspect > 1. ? param.aspect*base : base;
110 param.canvasY = param.aspect < 1. ? base/param.aspect : base;
111 } else {
112 param.canvasX = param.aspect > 1. ? param.aspect*1000 : T(1000);
113 param.canvasY = param.aspect < 1. ? 1000/param.aspect : T(1000);
114 }
115 // avoid exceeding jpeg pixel range
116 if (param.canvasX > 65500)
117 param.canvasX = 65499;
118 if (param.canvasY > 65500)
119 param.canvasY = 65499;
120 param.cbXscaling = param.canvasX/(double)1000;
122 writeHeatMapPlotFile(param, valueArea);
123 executeGnuplot(param);
124 }
125 }
126 return;
127}
128
129template <typename T>
131{
132
133
134 std::ofstream foutMatrix( param.matrixPath.c_str() );
135
136 int i[2] = {0,0};
137 for (i[1] = param.ny * param.zoomMin[1]; i[1] < param.ny * param.zoomMax[1]; i[1]++) {
138 for (i[0] = param.nx * param.zoomMin[0]; i[0] < param.nx * param.zoomMax[0]; i[0]++) {
139 T evaluated[1];
140 (*param.blockData)(evaluated,i);
141 foutMatrix << BaseType<T>(evaluated[0]) << " ";
142 if (i[0] == int(param.nx * param.zoomMax[0]) - 1) {
143 foutMatrix << "\n";
144 }
145 }
146 }
147 foutMatrix.close();
148
149
150 if (param.plot->writeCSV) {
151 std::ofstream foutCSV( param.csvPath.c_str() );
152 for (i[1] = 0; i[1] < param.ny; i[1]++) {
153 for (i[0] = 0; i[0] < param.nx; i[0]++) {
154 T evaluated[1];
155 Vector<T,3> physPoint = param.hyperPlane->getPhysR(i[0], i[1]);
156 (*param.blockData)(evaluated,i);
157 foutCSV << physPoint[0] << " " << physPoint[1] << " " << physPoint[2] << " " << evaluated[0] << "\n";
158 }
159 }
160 foutCSV.close();
161 }
162
163 return;
164}
165
166template <typename T>
167void writeHeatMapPlotFile(detailParam<T>& param, const std::vector<T>& valueArea)
168{
169 std::ofstream fout;
170
171 fout.open(param.plotFilePath.c_str(), std::ios::trunc);
172
173 fout << "if (strstrt(GPVAL_TERMINALS, 'jpeg') > 0) {";
174 fout << "set terminal jpeg " << "size " << param.canvasX << "," << param.canvasY << "font \",25\"" << "\n";
175 fout << "set output '"<< param.jpegPath << "'"<< "\n";
176 fout << "} else {";
177 fout << "set terminal png " << "size " << param.canvasX << "," << param.canvasY << "font \",25\"" << "\n";
178 fout << "set output '"<< param.pngPath << "'"<< "\n";
179 fout << "}" << "\n";
180 fout << "set pm3d map" << "\n";
181 fout << "unset key" << "\n";
182
183 if (param.plot->fullScreenPlot) {
184 fout << "unset xtics" << "\n";
185 fout << "unset ytics" << "\n";
186 fout << "unset border"<< "\n";
187 } else {
188 fout << "set xtics out" << "\n";
189 fout << "set ytics out" << "\n";
190 fout << "set xtics nomirror" << "\n";
191 fout << "set ytics nomirror" << "\n";
192 }
193
194
195
196 fout << "set pm3d interpolate 0,0" << "\n";
197
198 if (param.plot->fullScreenPlot) {
199 //plot without margins if fullscreenmode enabled
200 fout << "set lmargin at screen 0" << "\n";
201 fout << "set rmargin at screen 1" << "\n";
202 fout << "set tmargin at screen 0" << "\n";
203 fout << "set bmargin at screen 1" << "\n";
204 } else {
205 //set size ratio for non-fullscreenmode plots
206 fout << "set size ratio -1" << "\n";
207 fout << "set size 0.925,1.0" << "\n";
208 }
209
210 //enable colorbox
211 if( param.plot->fullScreenPlot == false || param.plot->activateFullScreenPlotColorBox == true ) {
212 fout << "set colorbox vertical user origin 0.85,0.1 size " << 0.025/BaseType<T>(param.cbXscaling) << " ,0.8" << "\n";
213 }
214
215 if ( util::nearZero(param.normal[0]) && util::nearZero(param.normal[1]) ) {
216 fout << "set xlabel \"x-axis in m \"" << "\n"
217 << "set ylabel \"y-axis in m \"" << "\n";
218 }
219 else if ( util::nearZero(param.normal[0]) && util::nearZero(param.normal[2]) ) {
220 fout << "set xlabel \"x-axis in m \"" << "\n"
221 << "set ylabel \"z-axis in m \"" << "\n";
222 param.origin[1] = param.origin[2];
223 }
224 else if ( util::nearZero(param.normal[1]) && util::nearZero(param.normal[2]) ) {
225 fout << "set xlabel \"y-axis in m \"" << "\n"
226 << "set ylabel \"z-axis in m \"" << "\n";
227 param.origin[0] = param.origin[1];
228 param.origin[1] = param.origin[2];
229 }
230 else {
231 fout << "set xlabel \"width in m \"" << "\n"
232 << "set ylabel \"height in m \"" << "\n";
233
234 }
235
236 if (param.plot->contourlevel > 0) {
237 fout << "set contour base" << "\n";
238 fout << "set cntrparam levels " << param.plot->contourlevel << "\n";
239 fout << "set cntrparam bspline" << "\n";
240 fout << "do for [i=1:"<< param.plot->contourlevel <<"] {" << "\n";
241 fout << "set linetype i lc rgb \"black\""<< "\n";
242 fout << "}" << "\n";
243 }
244
245 fout << "set cblabel offset 0.5 \"" <<param.quantityname;
246 if (param.quantityname == "l2(physVelocity)" || param.quantityname == "EuklidNorm(physVelocity)") {
247 fout << " in m/s\"" << "\n";
248 }
249 else if (param.quantityname == "physPressure") {
250 fout << " in Pa\"" << "\n";
251 }
252 else {
253 fout << "\"\n";
254 }
255
256 if (!util::nearZero(param.plot->maxValue - param.plot->minValue)) {
257 fout << "set cbrange [" << BaseType<T>(param.plot ->minValue) << ":" << BaseType<T>(param.plot->maxValue) <<"]" << "\n";
258 }
259
260 if (valueArea.empty()) {
261 fout << "set autoscale fix" << "\n";
262 }
263 else if (valueArea[0] < valueArea[1]) {
264 fout << "set cbrange [" << BaseType<T>(valueArea[0]) << ":" << BaseType<T>(valueArea[1]) << "]" << "\n";
265 }
266 else {
267 fout << "set cbrange [" << BaseType<T>(valueArea[1]) << ":" << BaseType<T>(valueArea[0]) << "]" << "\n";
268 }
269
270 if (param.plot->colour == "grey") {
271 fout << "set palette grey" << "\n";
272 }
273 else if (param.plot->colour == "pm3d") {
274
275 }
276 else if (param.plot->colour == "blackbody") {
277 fout << "set palette defined ( 0 \"black\", 1 \"red\", 2 \"yellow\")" << "\n";
278 }
279 else {
280 fout << "set palette defined ( 0 \"blue\", 1 \"green\", 2 \"yellow\", 3 \"orange\", 4 \"red\" )" << "\n";
281 }
282 fout << "splot '" << param.matrixPath << "' u ($1*" << BaseType<T>(param.spacing) << "+" << BaseType<T>(param.origin[0]) + int(param.nx * BaseType<T>(param.zoomMin[0]))*BaseType<T>(param.spacing) << "):"
283 << "($2*" << BaseType<T>(param.spacing) << "+" << BaseType<T>(param.origin[1]) + int(param.ny * BaseType<T>(param.zoomMin[1]))*BaseType<T>(param.spacing) <<"):3 matrix with pm3d" << "\n";
284 fout.close();
285 return;
286}
287
288template< typename T >
290{
291 if (! gnuplotInstalled()) {
292 std::cout << "We could not find a gnuplot distribution at your system." << std::endl;
293 std::cout << "We still write the data files s.t. you can plot the data yourself." << std::endl;
294 return;
295 }
296#ifndef WIN32
297 if (!system(nullptr)) {
298 exit (EXIT_FAILURE);
299 }
300 const std::string command = "gnuplot "+ param.plotFilePath +" > /dev/null &";
301 if ( system(command.c_str()) ) {
302 std::cout << "Error at GnuplotWriter" << std::endl;
303 }
304 return;
305#endif
306}
307
309{
310#ifdef WIN32
311 return false;
312#else
313 return (! system("which gnuplot >/dev/null 2>/dev/null"));
314#endif
315}
316
317} // namespace detail
318
319} // namespace heatmap
320
321} // namespace olb
322
323#endif
represents all functors that operate on a cuboid in general, mother class of BlockLatticeF,...
BlockReduction2D2D interpolates the data of a SuperF2D functor in a given resolution.
HyperplaneLattice3D< T > getPlaneDiscretizationIn3D() const
Returns embedding of the discretized plane in 3D space.
BlockReduction3D2D reduces the data of a SuperF3D functor to the intersection between a given hyperpl...
int getTargetDim() const
read only access to member variable _n
Definition genericF.hh:45
Parametrization of a hyperplane lattice.
Plain old scalar vector.
Definition vector.h:47
std::string getImageOutDir() const
Definition singleton.h:93
void writeHeatMapDataFile(detailParam< T > &param)
void genericHeatMapInterface(const HyperplaneLattice3D< T > &hyperPlane, BlockF2D< T > &blockData, int iT, const std::vector< T > &valueArea, const plotParam< T > &param)
void executeGnuplot(detailParam< T > &param)
void writeHeatMapPlotFile(detailParam< T > &param, const std::vector< T > &valueArea)
void write(BlockReduction3D2D< T > &blockReduction, int iT, const plotParam< T > param={}, const std::vector< T > &valueArea=std::vector< T >())
This function is used to plot heat maps as jpeg files.
MpiManager & mpi()
Directories & directories()
Definition singleton.h:150
bool nearZero(const ADf< T, DIM > &a)
Definition aDiff.h:1087
Top level namespace for all of OpenLB.
std::string createFileName(std::string name)
for .pvd masterFile
Definition fileName.hh:34
typename util::BaseTypeHelper< T >::type BaseType
Definition baseType.h:59
constexpr Vector< T, 3 > crossProduct3D(const ScalarVector< T, 3, IMPL > &a, const ScalarVector< T, 3, IMPL_ > &b) any_platform
Definition vector.h:224
Definition of singletons: global, publicly available information.
const HyperplaneLattice3D< T > * hyperPlane