MADNESS 0.10.1
funcplot.h
Go to the documentation of this file.
1/*
2 This file is part of MADNESS.
3
4 Copyright (C) 2007,2010 Oak Ridge National Laboratory
5
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19
20 For more information please contact:
21
22 Robert J. Harrison
23 Oak Ridge National Laboratory
24 One Bethel Valley Road
25 P.O. Box 2008, MS-6367
26
27 email: harrisonrj@ornl.gov
28 tel: 865-241-3937
29 fax: 865-572-0680
30*/
31
32#ifndef MADNESS_MRA_FUNCPLOT_H__INCLUDED
33#define MADNESS_MRA_FUNCPLOT_H__INCLUDED
34
35#include <madness/constants.h>
37/*!
38
39 \file mra/funcplot.h
40 \brief Defines/implements plotting interface for functions
41 \ingroup funcplot
42
43 @{
44 */
45
46namespace madness {
47
49 PlotParameters(World& world, const commandlineparser parser=commandlineparser(), const std::string tag="plot") : PlotParameters() {
50 read_input_and_commandline_options(world,parser,tag);
51 }
52
54
55 // initialize with: key, value, comment (optional), allowed values (optional)
56 initialize<double>("zoom",2,"zoom into the simulation cell");
57 initialize<long>("npoints",151,"number of plot points per dimension");
58 initialize<std::vector<double>>("origin",{},"origin of the plot");
59 initialize<std::vector<std::string>>("plane",{"x1","x2"},"plot plane: x1, x2, .., x6");
60 }
61
62 PlotParameters& set_zoom(const double z) {
63 set_user_defined_value("zoom",z);
64 return *this;
65 }
66 PlotParameters& set_npoints(const long n) {
67 set_user_defined_value("npoints",n);
68 return *this;
69 }
70 PlotParameters& set_plane(const std::vector<std::string> plane) {
72 return *this;
73 }
74 PlotParameters& set_origin(const std::vector<double> origin) {
76 return *this;
77 }
78
79
80 double zoom() const {return get<double>("zoom");}
81 long npoints() const {return get<long>("npoints");}
82
83 template<std::size_t NDIM>
85 auto origin_vec=get<std::vector<double>>("origin");
86 // fill in zeros if the default origin has fewer dimensions than the actual origin
87 int missing=NDIM-origin_vec.size();
88 for (auto i=0; i<missing; ++i) origin_vec.push_back(0.0);
89 Vector<double,NDIM> o(origin_vec);
90 return o;
91 }
92
93 std::vector<std::string> plane() const {return get<std::vector<std::string>>("plane");}
94
95
96 };
97 /// Writes an OpenDX format file with a cube/slice of points on a uniform grid
98
99 /// Collective operation but only process 0 writes the file. By convention OpenDX
100 /// files end in ".dx" but this choice is up to the user. The binary format is
101 /// more compact and vastly faster to both write and load but is not as portable.
102 ///
103 /// Now follow some brief tips about how to look at files inside OpenDX.
104 ///
105 /// To view a 1D function \c file-selector-->import-->plot-->image.
106 ///
107 /// To view a 2D function as a colored plane \c file-selector-->import-->autocolor-->image.
108 ///
109 /// To view a 2D function as a 3D surface \c file-selector-->import-->rubbersheet-->image.
110 ///
111 /// To view a 3D function as an isosurface \c file-selector-->import-->isosurface-->image.
112 ///
113 /// To select the real/imaginary/absolute value of a complex number insert a compute
114 /// element after the import.
115 template <typename T, std::size_t NDIM>
116 void plotdx(const Function<T,NDIM>& f,
117 const char* filename,
118 const Tensor<double>& cell = FunctionDefaults<NDIM>::get_cell(),
119 const std::vector<long>& npt = std::vector<long>(NDIM,201L),
120 bool binary=true);
121
122
123 /// Writes the header information of a VTK file for plotting in an external
124 /// post-processing package (such as Paraview)
125 //
126 /// @param world World communicator
127 /// @param filename String containing the filename to export to
128 /// @param plotlo Vector of double values indicating the minimum coordinate to plot to in each dimension
129 /// @param plothi Vector of double values indicating the maximum coordinate to plot to in each dimension
130 /// @param npt Vector of long integers indicating the number of points to plot in each dimension
131 /// @param binary (optional) Boolean indicating whether to print in binary
132
133 /// The VTK routines are also designed for SERIAL data, parallel coming...
134 ///
135 /// This header is templated by the dimension of the data.
136 ///
137 /// To plot with the plotvtk_* routines:
138 /// plotvtk_begin(...)
139 /// plotvtk_data(...)
140 /// plotvtk_data(...) ...
141 /// plotvtk_end(...)
142 ///
143 /// NOTE: Paraview expects the structured mesh points in a particular
144 /// order, which is why the LowDimIndexIterator is used...
145 template<std::size_t NDIM>
146 void plotvtk_begin(World &world, const char *filename,
147 const Vector<double, NDIM> &plotlo, const Vector<double, NDIM> &plothi,
148 const Vector<long, NDIM> &npt, bool binary = false) {
149
151 MADNESS_ASSERT(NDIM>=1 && NDIM<=3); // how do we plot data in more than 3-D?
152
153 Tensor<double> cell(NDIM, 2);
154 std::size_t i;
155 for(i = 0; i < NDIM; ++i) {
156 cell(i, 0) = plotlo[i];
157 cell(i, 1) = plothi[i];
158 }
159
160 FILE *f=0;
161 if(world.rank() == 0) {
162 f = fopen(filename, "w");
163 if(!f)
164 MADNESS_EXCEPTION("plotvtk: failed to open the plot file", 0);
165
166 fprintf(f, "<VTKFile type=\"StructuredGrid\" version=\"0.1\"" \
167 " byte_order=\"LittleEndian\" compressor=\"" \
168 "vtkZLibDataCompressor\">\n");
169 fprintf(f, " <StructuredGrid WholeExtent=\"");
170 for(i = 0; i < NDIM; ++i)
171 fprintf(f, "0 %ld ", npt[i]-1);
172 for(; i < 3; ++i)
173 fprintf(f, "0 0 ");
174 fprintf(f, "\">\n");
175 fprintf(f, " <Piece Extent=\"");
176 for(i = 0; i < NDIM; ++i)
177 fprintf(f, "0 %ld ", npt[i]-1);
178 for(; i < 3; ++i)
179 fprintf(f, "0 0 ");
180 fprintf(f, "\">\n");
181 fprintf(f, " <Points>\n");
182 fprintf(f, " <DataArray NumberOfComponents=\"3\" " \
183 "type=\"Float32\" format=\"ascii\">\n");
184
186 for(i = 0; i < NDIM; ++i) {
187 if(npt[i] == 1)
188 space[i] = 0.0;
189 else
190 space[i] = (cell(i, 1) - cell(i, 0)) / (npt[i] - 1);
191 }
192
193 // go through the grid
194 for(LowDimIndexIterator it(npt); it; ++it) {
195 for(i = 0; i < NDIM; ++i)
196 fprintf(f, "%f ", plotlo[i] + it[i]*space[i]);
197 for(; i < 3; ++i)
198 fprintf(f, "0.0 ");
199 fprintf(f, "\n");
200 }
201
202 fprintf(f, " </DataArray>\n");
203 fprintf(f, " </Points>\n");
204 fprintf(f, " <PointData>\n");
205 fclose(f);
206 }
207 world.gop.fence();
208 }
209
210 /// Generic VTK data writer. Specific type instances of this function are defined for
211 /// both real and complex valued functions.
212 //
213 /// @param function Function (real or complex) that we wish to export the data of
214 /// @param fieldname A string containing the name we wish to refer to this field as in the exported data
215 /// @param world World communicator
216 /// @param filename String containing the filename to export to
217 /// @param plotlo Vector of double values indicating the minimum coordinate to plot to in each dimension
218 /// @param plothi Vector of double values indicating the maximum coordinate to plot to in each dimension
219 /// @param npt Vector of long integers indicating the number of points to plot in each dimension
220 /// @param binary (optional) Boolean indicating whether to print in binary
221
222 /// This templated function won't do anything except print a warning
223 /// message. Specialized versions of this function should be used.
224 template<typename T, std::size_t NDIM>
225 void plotvtk_data(const T &function, const char *fieldname, World &world,
226 const char *filename, const Vector<double, NDIM> &plotlo,
227 const Vector<double, NDIM> &plothi, const Vector<long, NDIM> &npt,
228 bool binary = false) {
229
230 MADNESS_EXCEPTION("plotvtk only supports madness::functions", 0);
231 }
232
233 /// VTK data writer for real-valued (not complex) madness::functions.
234
235 /// Set plot_refine=true to get a plot of the refinement levels of
236 /// the given function.
237 template<typename T, std::size_t NDIM>
238 void plotvtk_data(const Function<T, NDIM> &function, const char *fieldname,
239 World &world, const char *filename, const Vector<double, NDIM> &plotlo,
240 const Vector<double, NDIM> &plothi, const Vector<long, NDIM> &npt,
241 bool binary = false, bool plot_refine = false) {
242
244 MADNESS_ASSERT(NDIM>=1 && NDIM<=3); // no plotting high-D functions, yet...
245
246 Tensor<double> cell(NDIM, 2);
247 std::size_t i;
248 for(i = 0; i < NDIM; ++i) {
249 cell(i, 0) = plotlo[i];
250 cell(i, 1) = plothi[i];
251 }
252 std::vector<long> numpt(NDIM);
253 for(i = 0; i < NDIM; ++i)
254 numpt[i] = npt[i];
255
256 world.gop.barrier();
257
258 function.verify();
259 FILE *f = 0;
260 if(world.rank() == 0) {
261 f = fopen(filename, "a");
262 if(!f)
263 MADNESS_EXCEPTION("plotvtk: failed to open the plot file", 0);
264
265 fprintf(f, " <DataArray Name=\"%s\" format=\"ascii\" " \
266 "type=\"Float32\" NumberOfComponents=\"1\">\n", fieldname);
267 }
268
269 world.gop.fence();
270 Tensor<T> tmpr = function.eval_cube(cell, numpt, plot_refine);
271 world.gop.fence();
272
273 if(world.rank() == 0) {
274 for(LowDimIndexIterator it(numpt); it; ++it) {
275 fprintf(f, "%.6e\n", tmpr(*it));
276 }
277 fprintf(f, " </DataArray>\n");
278 fclose(f);
279 }
280 world.gop.fence();
281 }
282
283 /// VTK data writer for complex-valued madness::functions.
284
285 /// The complex-value is written as two reals (a vector from VTK's
286 /// perspective. The first (X) component is the real part and the second
287 /// (Y) component is the imaginary part.
288 /// Set plot_refine=true to get a plot of the refinement levels of
289 /// the given function.
290 template<typename T, std::size_t NDIM>
291 void plotvtk_data(const Function<std::complex<T>, NDIM> &function,
292 const char *fieldname, World &world, const char *filename,
293 const Vector<double, NDIM> &plotlo, const Vector<double, NDIM> &plothi,
294 const Vector<long, NDIM> &npt, bool binary = false,
295 bool plot_refine = false) {
296
297 // this is the same as plotvtk_data for real functions, except the
298 // real and imaginary parts are printed on the same line (needed
299 // to change NumberOfComponents in the XML tag)
300
302 MADNESS_ASSERT(NDIM>=1 && NDIM<=3); // no plotting high-D functions, yet...
303
304 Tensor<double> cell(NDIM, 2);
305 std::size_t i;
306 for(i = 0; i < NDIM; ++i) {
307 cell(i, 0) = plotlo[i];
308 cell(i, 1) = plothi[i];
309 }
310 std::vector<long> numpt(NDIM);
311 for(i = 0; i < NDIM; ++i)
312 numpt[i] = npt[i];
313
314 world.gop.barrier();
315
316 function.verify();
317 FILE *f = 0;
318 if(world.rank() == 0) {
319 f = fopen(filename, "a");
320 if(!f)
321 MADNESS_EXCEPTION("plotvtk: failed to open the plot file", 0);
322
323 fprintf(f, " <DataArray Name=\"%s\" format=\"ascii\" " \
324 "type=\"Float32\" NumberOfComponents=\"2\">\n", fieldname);
325 }
326
327 world.gop.fence();
328 Tensor<std::complex<T> > tmpr = function.eval_cube(cell, numpt,
329 plot_refine);
330 world.gop.fence();
331
332 if(world.rank() == 0) {
333 for(LowDimIndexIterator it(numpt); it; ++it) {
334 fprintf(f, "%.6e %.6e\n", real(tmpr(*it)), imag(tmpr(*it)));
335 }
336 fprintf(f, " </DataArray>\n");
337 fclose(f);
338 }
339 world.gop.fence();
340 }
341
342 /// Writes the footer information of a VTK file for plotting in an external
343 /// post-processing package (such as Paraview)
344 //
345 /// @param world World communicator
346 /// @param filename Name of VTK file
347 /// @param binary (Optional) Boolean indicating whether to print in binary
348 template<std::size_t NDIM>
349 void plotvtk_end(World &world, const char *filename, bool binary = false) {
351 MADNESS_ASSERT(NDIM>=1 && NDIM<=3);
352
353 FILE *f = 0;
354 if(world.rank() == 0) {
355 f = fopen(filename, "a");
356 if(!f)
357 MADNESS_EXCEPTION("plotvtk: failed to open the plot file", 0);
358
359 fprintf(f, " </PointData>\n");
360 fprintf(f, " <CellData>\n");
361 fprintf(f, " </CellData>\n");
362 fprintf(f, " </Piece>\n");
363 fprintf(f, " </StructuredGrid>\n");
364 fprintf(f, "</VTKFile>\n");
365 fclose(f);
366 }
367 world.gop.fence();
368 }
369
370 namespace detail {
371 inline unsigned short htons_x(unsigned short a) {
372 return (a>>8) | (a<<8);
373 }
374 }
375
376 /// Writes a Povray DF3 format file with a cube of points on a uniform grid
377
378 /// Collective operation but only process 0 writes the file. By convention Povray
379 /// files end in ".df3" but this choice is up to the user. The dynamic range of
380 /// function values is mapped onto [0,1] and values stored in 16-bit fixed precision.
381 template <typename T>
383 const char* filename,
385 const std::vector<long>& npt = std::vector<long>(3,201L))
386 {
387 using detail::htons_x;
388
389 MADNESS_ASSERT(npt.size() == 3);
390 unsigned short dims[3] = {htons_x(npt[0]),htons_x(npt[1]),htons_x(npt[2])};
391
392 World& world = const_cast< Function<T,3>& >(function).world();
393 FILE *f=0;
394 if (world.rank() == 0) {
395 f = fopen(filename, "w");
396 if (!f) MADNESS_EXCEPTION("plotdx: failed to open the plot file", 0);
397 fwrite((void*) dims, sizeof(short), 3, f);
398 }
399 Tensor<T> r = function.eval_cube(cell, npt);
400 if (world.rank() == 0) {
401 double rmax = r.max();
402 double rmin = r.min();
403 double rrange = rmax + rmin;
404 double rmean = rrange*0.5;
405 double fac = 65535.0/rrange;
406
407 printf("plot_povray: %s: min=%.2e(0.0) mean=%.2e(0.5) max=%.2e(1.0) range=%.2e\n",
408 filename,rmin,rmean,rmax,rrange);
409
410 std::vector<unsigned short> d(npt[0]);
411 for (unsigned int i2=0; i2<npt[2]; ++i2) {
412 for (unsigned int i1=0; i1<npt[1]; ++i1) {
413 for (unsigned int i0=0; i0<npt[0]; ++i0) {
414 d[i0] = (unsigned short)(htons_x((unsigned short)(fac*(r(i0,i1,i2) - rmin))));
415 //printf("%d\n",htons_x(d[i0]));
416 }
417 fwrite((void*) &d[0], sizeof(short), npt[0], f);
418 }
419 }
420
421 fclose(f);
422 }
423 }
424
425 static inline void plot_line_print_value(FILE* f, double_complex v) {
426 fprintf(f, " %.14e %.14e ", real(v), imag(v));
427 }
428
429 static inline void plot_line_print_value(FILE* f, double v) {
430 fprintf(f, " %.14e", v);
431 }
432
433
434 /// Generates ASCII file tabulating f(r) at npoints along line r=lo,...,hi
435
436 /// The ordinate is distance from lo
437 template <typename opT, std::size_t NDIM>
438 void plot_line(World& world, const char* filename, int npt, const Vector<double,NDIM>& lo,
439 const Vector<double,NDIM>& hi, const opT& op) {
441 coordT h = (hi - lo)*(1.0/(npt-1));
442
443 double sum = 0.0;
444 for (std::size_t i=0; i<NDIM; ++i) sum += h[i]*h[i];
445 sum = sqrt(sum);
446
447 if (world.rank() == 0) {
448 FILE* file = fopen(filename,"w");
449 if(!file)
450 MADNESS_EXCEPTION("plot_line: failed to open the plot file", 0);
451 for (int i=0; i<npt; ++i) {
452 coordT r = lo + h*double(i);
453 fprintf(file, "%.14e ", i*sum);
454 plot_line_print_value(file, op(r));
455 fprintf(file,"\n");
456 }
457 fclose(file);
458 }
459 world.gop.fence();
460 }
461
462
463 /// Generates ASCII file tabulating f(r) at npoints along line r=lo,...,hi
464
465 /// The ordinate is distance from lo
466 template <typename T, std::size_t NDIM>
467 void plot_line(const char* filename, int npt, const Vector<double,NDIM>& lo, const Vector<double,NDIM>& hi,
468 const Function<T,NDIM>& f) {
470 coordT h = (hi - lo)*(1.0/(npt-1));
471
472 double sum = 0.0;
473 for (std::size_t i=0; i<NDIM; ++i) sum += h[i]*h[i];
474 sum = sqrt(sum);
475
476 World& world = f.world();
477 f.reconstruct();
478 if (world.rank() == 0) {
479 FILE* file = fopen(filename,"w");
480 if(!file)
481 MADNESS_EXCEPTION("plot_line: failed to open the plot file", 0);
482 for (int i=0; i<npt; ++i) {
483 coordT r = lo + h*double(i);
484 fprintf(file, "%.14e ", i*sum);
485 plot_line_print_value(file, f.eval(r));
486 fprintf(file,"\n");
487 }
488 fclose(file);
489 }
490 world.gop.fence();
491 }
492
493 /// Generates ASCII file tabulating f(r) and g(r) at npoints along line r=lo,...,hi
494
495 /// The ordinate is distance from lo
496 template <typename T, typename U, std::size_t NDIM>
497 void plot_line(const char* filename, int npt, const Vector<double,NDIM>& lo, const Vector<double,NDIM>& hi,
498 const Function<T,NDIM>& f, const Function<U,NDIM>& g) {
500 coordT h = (hi - lo)*(1.0/(npt-1));
501
502 double sum = 0.0;
503 for (std::size_t i=0; i<NDIM; ++i) sum += h[i]*h[i];
504 sum = sqrt(sum);
505
506 World& world = f.world();
507 f.reconstruct();
508 g.reconstruct();
509 if (world.rank() == 0) {
510 FILE* file = fopen(filename,"w");
511 if(!file)
512 MADNESS_EXCEPTION("plot_line: failed to open the plot file", 0);
513 for (int i=0; i<npt; ++i) {
514 coordT r = lo + h*double(i);
515 fprintf(file, "%.14e ", i*sum);
516 plot_line_print_value(file, f.eval(r));
517 plot_line_print_value(file, g.eval(r));
518 fprintf(file,"\n");
519 }
520 fclose(file);
521 }
522 world.gop.fence();
523 }
524
525
526 /// Generates ASCII file tabulating f(r), g(r), and a(r) at npoints along line r=lo,...,hi
527
528 /// The ordinate is distance from lo
529 template <typename T, typename U, typename V, std::size_t NDIM>
530 void plot_line(const char* filename, int npt, const Vector<double,NDIM>& lo, const Vector<double,NDIM>& hi,
531 const Function<T,NDIM>& f, const Function<U,NDIM>& g, const Function<V,NDIM>& a) {
533 coordT h = (hi - lo)*(1.0/(npt-1));
534
535 double sum = 0.0;
536 for (std::size_t i=0; i<NDIM; ++i) sum += h[i]*h[i];
537 sum = sqrt(sum);
538
539 World& world = f.world();
540 f.reconstruct();
541 g.reconstruct();
542 a.reconstruct();
543 if (world.rank() == 0) {
544 FILE* file = fopen(filename,"w");
545 if(!file)
546 MADNESS_EXCEPTION("plot_line: failed to open the plot file", 0);
547 for (int i=0; i<npt; ++i) {
548 coordT r = lo + h*double(i);
549 fprintf(file, "%.14e ", i*sum);
550 plot_line_print_value(file, f.eval(r));
551 plot_line_print_value(file, g.eval(r));
552 plot_line_print_value(file, a.eval(r));
553 fprintf(file,"\n");
554 }
555 fclose(file);
556 }
557 world.gop.fence();
558 }
559
560 /// Generates ASCII file tabulating f(r), g(r), a(r), b(r) at npoints along line r=lo,...,hi
561
562 /// The ordinate is distance from lo
563 template <typename T, typename U, typename V, typename W, std::size_t NDIM>
564 void plot_line(const char* filename, int npt, const Vector<double,NDIM>& lo, const Vector<double,NDIM>& hi,
565 const Function<T,NDIM>& f, const Function<U,NDIM>& g, const Function<V,NDIM>& a, const Function<W,NDIM>& b) {
567 coordT h = (hi - lo)*(1.0/(npt-1));
568
569 double sum = 0.0;
570 for (std::size_t i=0; i<NDIM; ++i) sum += h[i]*h[i];
571 sum = sqrt(sum);
572
573 World& world = f.world();
574 f.reconstruct();
575 g.reconstruct();
576 a.reconstruct();
577 b.reconstruct();
578 if (world.rank() == 0) {
579 FILE* file = fopen(filename,"w");
580 for (int i=0; i<npt; ++i) {
581 coordT r = lo + h*double(i);
582 fprintf(file, "%.14e ", i*sum);
583 plot_line_print_value(file, f.eval(r));
584 plot_line_print_value(file, g.eval(r));
585 plot_line_print_value(file, a.eval(r));
586 plot_line_print_value(file, b.eval(r));
587 fprintf(file,"\n");
588 }
589 fclose(file);
590 }
591 world.gop.fence();
592 }
593 /// The ordinate is distance from lo
594 template <typename T, std::size_t NDIM>
595 void plot_line(const char* filename, int npt, const Vector<double,NDIM>& lo, const Vector<double,NDIM>& hi,
596 const std::vector<Function<T,NDIM>>& vf) {
598 coordT h = (hi - lo)*(1.0/(npt-1));
599 double sum = 0.0;
600 for (std::size_t i=0; i<NDIM; ++i) sum += h[i]*h[i];
601 sum = sqrt(sum);
602 World& world = vf[0].world();// get world from first function
603 // reconstruct each function in vf
604 std::for_each(vf.begin(), vf.end(), [](const Function<T,NDIM>& f){f.reconstruct();});
605 if (world.rank() == 0) {
606 FILE* file = fopen(filename,"w");
607 if(!file)
608 MADNESS_EXCEPTION("plot_line: failed to open the plot file", 0);
609 for (int i=0; i<npt; ++i) {
610 coordT r = lo + h*double(i);
611 fprintf(file, "%.14e ", i*sum);
612 std::for_each(vf.begin(), vf.end(), [&](const Function<T,NDIM>& f){ plot_line_print_value(file, f.eval(r));});
613 fprintf(file,"\n");
614 }
615 fclose(file);
616 }
617 world.gop.fence();
618 }
619
620 template<size_t NDIM>
622 const std::string name) {
623 typedef std::vector<Function<double,NDIM> > vecfuncT;
625 }
626
627 template<size_t NDIM>
628 void plot_plane(World& world, const Function<double,NDIM>& function1,
629 const Function<double,NDIM>& function2,
630 const std::string name) {
631 typedef std::vector<Function<double,NDIM> > vecfuncT;
632 vecfuncT vf(2);
633 vf[0]=function1;
634 vf[1]=function2;
635 plot_plane(world,vf,name);
636 }
637
638 template<size_t NDIM>
639 void plot_plane(World& world, const Function<double,NDIM>& function1,
640 const Function<double,NDIM>& function2,const Function<double,NDIM>& function3,
641 const std::string name) {
642 typedef std::vector<Function<double,NDIM> > vecfuncT;
643 vecfuncT vf(3);
644 vf[0]=function1;
645 vf[1]=function2;
646 vf[2]=function3;
647 plot_plane(world,vf,name);
648 }
649
650
651 /// plot a 2-d slice of a given function and the according MRA structure
652 /// FIXME: doesn't work for more than 1 rank
653
654 /// the plotting parameters are taken from the input file "input" and its
655 /// data group "plot", e.g. plotting the xy plane around (0,0,0.7):
656 /// plot
657 /// plane x1 x2
658 /// zoom 2.0
659 /// npoints 100
660 /// origin 0.0 0.0 0.7
661 /// end
662 /// @param[in] world the world
663 /// @param[in] vfunction the function to plot
664 /// @param[in] name the output name
665 template<size_t NDIM>
666 void plot_plane(World& world, const std::vector<Function<double,NDIM> >& vfunction,
667 const std::string name, const std::string inputfile="input") {
668
669 if (world.size()>1) return;
670 // determine the ploting plane
671 std::string c1="x1", c2="x2";
672
673 // zoom factor
674 double zoom=1.0;
675
676 // output type: mathematica or gnuplot
677 std::string output_type="gnuplot";
678
679 // number of points in each direction
680 int npoints=200;
681
682 // the coordinates to be plotted
683 Vector<double,NDIM> coord(0.0);
684 Vector<double,NDIM> origin(0.0);
685
686 try {
687 std::ifstream f(inputfile);
688 position_stream_to_word(f, "plot",'#',true,true);
689 std::string s;
690 while (f >> s) {
691 if (s == "end") {
692 break;
693 } else if (s == "plane") {
694 f >> c1 >> c2;
695 } else if (s == "zoom") {
696 f >> zoom;
697 } else if (s == "output") {
698 f >> output_type;
699 } else if (s == "points") {
700 f >> npoints;
701 } else if (s == "origin") {
702 for (std::size_t i=0; i<NDIM; ++i) f >> origin[i];
703 }
704 }
705 } catch (...) {
706 print("can't locate plot in file="+inputfile+" -- using default values");
707 }
708 double scale=1.0/zoom;
709 coord=origin;
710
711 // convert human to mad form
712 size_t cc1=0, cc2=1;
713 if (c1=="x1") cc1=0;
714 if (c1=="x2") cc1=1;
715 if (c1=="x3") cc1=2;
716 if (c1=="x4") cc1=3;
717 if (c1=="x5") cc1=4;
718 if (c1=="x6") cc1=5;
719 if (c2=="x1") cc2=0;
720 if (c2=="x2") cc2=1;
721 if (c2=="x3") cc2=2;
722 if (c2=="x4") cc2=3;
723 if (c2=="x5") cc2=4;
724 if (c2=="x6") cc2=5;
725
726 MADNESS_ASSERT(cc1<NDIM);
727 MADNESS_ASSERT(cc2<NDIM);
728 // output file name for the gnuplot data
729 std::string filename="plane_"+c1+c2+"_"+name;
730 // assume a cubic cell
732 lo=lo*scale;
733
734 const double stepsize=FunctionDefaults<NDIM>::get_cell_width()[0]*scale/npoints;
735
736 if(world.rank() == 0) {
737
738 // plot 3d plot
739 FILE *f = 0;
740 f=fopen(filename.c_str(), "w");
741 if(!f) MADNESS_EXCEPTION("plot_along: failed to open the plot file", 0);
742
743 for (int i0=0; i0<npoints; i0++) {
744 for (int i1=0; i1<npoints; i1++) {
745 // plot plane
746 coord[cc1]=lo+origin[cc1]+i0*stepsize;
747 coord[cc2]=lo+origin[cc2]+i1*stepsize;
748
749 // other electron
750// fprintf(f,"%12.6f %12.6f %12.20f\n",coord[cc1],coord[cc2],
751// function(coord));
752 fprintf(f,"%12.6f %12.6f",coord[cc1],coord[cc2]);
753 for (std::size_t ivec=0; ivec<vfunction.size(); ++ivec)
754 fprintf(f," %12.20f",vfunction[ivec](coord));
755 fprintf(f,"\n");
756
757 }
758 // additional blank line between blocks for gnuplot
759 if (output_type=="gnuplot") fprintf(f,"\n");
760 }
761 fclose(f);
762
763 }
764
765// // plot mra structure
766// filename="mra_structure_"+c1+c2+"_"+name;
767// function.get_impl()->print_plane(filename.c_str(),cc1,cc2,coord);
768 }
769
770
771/// plot a 2-d slice of a given function and the according MRA structure
772
773/// the plotting parameters are taken from the input file "input" and its
774/// data group "plot", e.g. plotting the xy plane around (0,0,0.7):
775/// plot
776/// plane x1 x2
777/// zoom 2.0
778/// npoints 100
779/// origin 0.0 0.0 0.7
780/// end
781/// @param[in] world the world
782/// @param[in] vfunction the function to plot
783/// @param[in] name the output name
784template<size_t NDIM>
785void plot_plane(World& world, const std::vector<Function<double,NDIM> >& vfunction,
786 const std::string name, const PlotParameters param) {
787
788 if (world.size()>1) return;
789
790 auto plane=param.plane();
791 std::string c1=plane[0];
792 std::string c2=plane[1];
793 auto npoints=param.npoints();
794 auto origin=param.origin<NDIM>();
795 auto coord=param.origin<NDIM>();
796 double scale=1.0/param.zoom();
797 std::string output_type="gnuplot";
798
799 auto plane2dim = [](std::string c) {
800 if (c=="x1") return 0;
801 else if (c=="x2") return 1;
802 else if (c=="x3") return 2;
803 else if (c=="x4") return 3;
804 else if (c=="x5") return 4;
805 else if (c=="x6") return 5;
806 else return -1;
807 };
808
809 // convert human to mad form
810 std::size_t cc1=plane2dim(c1);
811 std::size_t cc2=plane2dim(c2);
812
813 MADNESS_ASSERT(cc1>=0 && cc1<NDIM);
814 MADNESS_ASSERT(cc2>=0 && cc2<NDIM);
815
816 // output file name for the gnuplot data
817 std::string filename="plane_"+c1+c2+"_"+name;
818 // assume a cubic cell
820 lo=lo*scale;
821
822 const double stepsize=FunctionDefaults<NDIM>::get_cell_width()[0]*scale/npoints;
823
824 if(world.rank() == 0) {
825
826 // plot 3d plot
827 FILE *f = 0;
828 f=fopen(filename.c_str(), "w");
829 if(!f) MADNESS_EXCEPTION("plot_along: failed to open the plot file", 0);
830
831 for (int i0=0; i0<npoints; i0++) {
832 for (int i1=0; i1<npoints; i1++) {
833 // plot plane
834 coord[cc1]=lo+origin[cc1]+i0*stepsize;
835 coord[cc2]=lo+origin[cc2]+i1*stepsize;
836
837 fprintf(f,"%12.6f %12.6f",coord[cc1],coord[cc2]);
838 for (std::size_t ivec=0; ivec<vfunction.size(); ++ivec)
839 fprintf(f," %12.20f",vfunction[ivec](coord));
840 fprintf(f,"\n");
841
842 }
843 // additional blank line between blocks for gnuplot
844 if (output_type=="gnuplot") fprintf(f,"\n");
845 }
846 fclose(f);
847
848 }
849}
850
851
852 template<size_t NDIM, typename opT>
853 void plot_plane(World& world, const opT& op, const std::string name, const PlotParameters param) {
854
855 if (world.size()>1) return;
856
857 auto plane=param.plane();
858 std::string c1=plane[0];
859 std::string c2=plane[1];
860 auto npoints=param.npoints();
861 auto origin=param. template origin<NDIM>();
862 auto coord=param. template origin<NDIM>();
863 double scale=1.0/param.zoom();
864 std::string output_type="gnuplot";
865
866 auto plane2dim = [](std::string c) {
867 if (c=="x1") return 0;
868 else if (c=="x2") return 1;
869 else if (c=="x3") return 2;
870 else if (c=="x4") return 3;
871 else if (c=="x5") return 4;
872 else if (c=="x6") return 5;
873 else return -1;
874 };
875
876 // convert human to mad form
877 std::size_t cc1=plane2dim(c1);
878 std::size_t cc2=plane2dim(c2);
879
880 MADNESS_ASSERT(cc1>=0 && cc1<NDIM);
881 MADNESS_ASSERT(cc2>=0 && cc2<NDIM);
882
883 // output file name for the gnuplot data
884 std::string filename="plane_"+c1+c2+"_"+name;
885 // assume a cubic cell
887 lo=lo*scale;
888
889 const double stepsize=FunctionDefaults<NDIM>::get_cell_width()[0]*scale/npoints;
890
891 if(world.rank() == 0) {
892
893 // plot 3d plot
894 FILE *f = 0;
895 f=fopen(filename.c_str(), "w");
896 if(!f) MADNESS_EXCEPTION("plot_along: failed to open the plot file", 0);
897
898 for (int i0=0; i0<npoints; i0++) {
899 for (int i1=0; i1<npoints; i1++) {
900 // plot plane
901 coord[cc1]=lo+origin[cc1]+i0*stepsize;
902 coord[cc2]=lo+origin[cc2]+i1*stepsize;
903
904 // other electron
905 // fprintf(f,"%12.6f %12.6f %12.20f\n",coord[cc1],coord[cc2],
906 // function(coord));
907 fprintf(f,"%12.6f %12.6f",coord[cc1],coord[cc2]);
908 fprintf(f," %12.20f\n",op(coord));
909
910 }
911 // additional blank line between blocks for gnuplot
912 if (output_type=="gnuplot") fprintf(f,"\n");
913 }
914 fclose(f);
915
916 }
917 }
918
919
920
921 template<size_t NDIM>
922 typename std::enable_if<NDIM==3,void>::type
924 std::vector<std::string> molecular_info=std::vector<std::string>(), int npoints=100, double zoom=1.0) {
925
926 if (world.size()>1) return;
927
928 // dummy atom in the center
929 if (molecular_info.size()==0)
930 molecular_info=std::vector<std::string>(1,"0 0 0.0 0.0 0.0\n");
931
932 // the coordinates to be plotted
933 Vector<double,NDIM> origin(0.0);
934
935 // number of points in each direction
936 std::vector<int> npt(3,npoints);
937
939 cell.scale(1.0/zoom);
940 double xlen=cell(0,1)-cell(0,0);
941 double ylen=cell(1,1)-cell(1,0);
942 double zlen=cell(2,1)-cell(2,0);
943
944 // plot file
945 FILE *file = 0;
946 file=fopen(filename.c_str(), "w");
947 if(!file) MADNESS_EXCEPTION("plot_along: failed to open the plot file", 0);
948
949
950 // print header
951 fprintf(file,"cube file from MADNESS\n");
952 fprintf(file,"comment line\n");
953
954 // print the number of atoms if a calculation was provided
955 fprintf(file,"%d %12.8f %12.8f %12.8f \n",int(molecular_info.size()),
956 cell(0,0),cell(1,0),cell(2,0));
957
958 // grid spacing for each dimension such that the cell edges are plotted
959 const auto xdelta = xlen/(npt[0]-1);
960 const auto ydelta = ylen/(npt[1]-1);
961 const auto zdelta = zlen/(npt[2]-1);
962
963 // print the cell constants
964 fprintf(file,"%d %12.6f %12.6f %12.6f\n",npt[0],xdelta,0.0,0.0);
965 fprintf(file,"%d %12.6f %12.6f %12.6f\n",npt[1],0.0,ydelta,0.0);
966 fprintf(file,"%d %12.6f %12.6f %12.6f\n",npt[2],0.0,0.0,zdelta);
967
968 // print the molecule
969 for (const std::string& s : molecular_info) fprintf(file,"%s",s.c_str());
970
971
972 Tensor<double> grid(npt[0], npt[1], npt[2]);
973 long count_per_line = 0;
974 for (int i = 0; i < npt[0]; ++i) {
975 for (int j = 0; j < npt[1]; ++j) {
976 for (int k = 0; k < npt[2]; ++k) {
977 double x = cell(0, 0) + origin[0] + xdelta * i;
978 double y = cell(1, 0) + origin[1] + ydelta * j;
979 double z = cell(2, 0) + origin[2] + zdelta * k;
980 // the original format has up to 6 entries per line: https://paulbourke.net/dataformats/cube/
981 // stick with this, even though many codes can read an arbitrary number of entries per line
982 if (count_per_line < 6) {
983 fprintf(file, "%12.5e ", f(x, y, z));
984 count_per_line++;
985 } else {
986 fprintf(file, "%12.5e\n", f(x, y, z));
987 count_per_line = 0;
988 }
989 }
990 }
991 }
992 fprintf(file, "\n");
993 fclose(file);
994 }
995
996 template<size_t NDIM>
997 typename std::enable_if<NDIM==3,void>::type
999
1000 if (world.size() > 1)
1001 return;
1002
1004 std::ofstream os(filename.c_str());
1005
1006 os << "{";
1007 os << "\"cell\":[";
1008 for (int xyz = 0; xyz != 3; ++xyz) {
1009 os << "[" << cell(xyz, 0) << "," << cell(xyz, 1) << "]";
1010 if (xyz != 2)
1011 os << ",";
1012 }
1013 os << "],";
1014
1015 os << "\"tree\":{";
1016 f.print_tree_json(os);
1017 os << "}}";
1018 }
1019
1020 /// convenience to get plot_plane and plot_cubefile
1021 template<size_t NDIM>
1022 void plot(const std::vector<Function<double,NDIM> >& vf, const std::string& name, const std::vector<std::string>& header){
1023 if(vf.empty()) return;
1024 World& world=vf.front().world();
1025 for(size_t i=0;i<vf.size();++i){
1026 const std::string namei=name+"_"+std::to_string(i);
1027// vf[i].print_size("plot:"+namei);
1028 plot_plane<NDIM>(world,vf[i],namei);
1029// plot_cubefile<NDIM>(world,vf[i],namei+".cube",header);
1030 }
1031 }
1032
1033 template<typename T>
1034 static std::string stringify(T arg) {
1035 std::ostringstream o;
1036 if (!(o << arg))
1037 MADNESS_EXCEPTION("stringify<T> failed",1);
1038 return o.str();
1039 }
1040
1041
1044
1045 // plot along this trajectory
1046 template<size_t NDIM>
1047 struct trajectory {
1048
1050
1051 double lo;
1052 double hi;
1053 double radius;
1054 long npt;
1055 coordT start, end;
1057 coordT (*curve)(const coordT& lo, const coordT& hi, double radius, coord_3d el2, long npt, long ipt);
1058
1059 /// some tools for plotting MRA ranks of low order tensors
1060
1061 // return a hue number [0,0.7] according to the rank in relation to maxrank,
1062 static double hueCode(const int rank) {
1063 const double maxrank=40.0;
1064 double hue=0.7-(0.7/maxrank)*(rank);
1065 return std::max(0.0,hue);
1066 }
1067
1068
1069 // print a dot of hue color at (x,y) in file f
1070 static void print_psdot(FILE *f, double x, double y, double color) {
1071 fprintf(f,"\\newhsbcolor{mycolor}{%8.4f 1.0 0.7}\n",color);
1072 fprintf(f,"\\psdot[linecolor=mycolor](%12.8f,%12.8f)\n",x,y);
1073 }
1074
1075
1076 static coord_3d circle2(double lo, double hi, double radius, coord_3d el2, long npt, long ipt) {
1077 double stepsize=constants::pi * 2.0 / npt;
1078 double phi=ipt*stepsize;
1079
1080 // in the xz plane
1081 coord_3d coord(0.0);
1082 coord[0]=radius * sin(phi);
1083 coord[1]=radius * cos(phi);
1084 return coord;
1085
1086 }
1087
1088 static coord_6d circle_6d(const coord_6d& lo, const coord_6d& hi, double radius, coord_3d el2, long npt, long ipt) {
1089 double stepsize=constants::pi * 2.0 / npt;
1090
1091 // start at phi=1.0
1092 double phi=1.0+constants::pi+ipt*stepsize;
1093
1094 // in the xz plane
1095 coord_6d coord(0.0);
1096 coord[0]=radius * sin(phi);
1097 coord[1]=radius * cos(phi);
1098 coord[2]=0.0;
1099 coord[3]=el2[0];
1100 coord[4]=el2[1];
1101 coord[5]=el2[2];
1102
1103 return coord;
1104
1105 }
1106
1107
1108 // typedef Vector<double,NDIM> (trajectory::circle_6d)(double lo, double hi, double radius, long npt, long ipt) const;
1109
1111// // ctor for a straight line thru the origin
1112// trajectory(double lo, double hi, long npt) : lo(lo), hi(hi), npt(npt), curve(line) {
1113// }
1114
1115 // ctor for circle
1116 trajectory(double radius, long npt) : radius(radius), npt(npt), curve(this->circle2) {
1117 }
1118
1119 // ctor for circle with electron 2 fixed at coord_3d
1120 trajectory(double radius, coord_3d el2, long npt) : radius(radius), npt(npt), el2(el2), curve(this->circle_6d) {
1121 }
1122
1123
1124 static Vector<double, NDIM> line_internal(const coordT& lo, const coordT& hi, double radius, coord_3d el2, long npt, long ipt) {
1125 const coordT step=(hi-lo)*(1.0/npt);
1126 coordT coord=lo+step*ipt;
1127 return coord;
1128 }
1129
1130
1131 /// constructor for a line
1132// static trajectory line(const Vector<double,NDIM>& lo, const Vector<double,NDIM>& hi, const long npt) {
1133 static trajectory line2(const coordT start, const coordT end, const long npt) {
1134 trajectory<NDIM> traj;
1135 traj.start=start;
1136 traj.end=end;
1137 traj.npt=npt;
1139 return traj;
1140 }
1141
1142 /// EZ ctor for a line a direction xyz={0,1,2,..,NDIM-1} thru the origin
1143 static trajectory line_xyz(const int xyz, const long npt) {
1145 coordT lo(0.0), hi(0.0);
1146 lo[xyz]=-L/2;
1147 hi[xyz]=L/2;
1148 return trajectory<NDIM>::line2(lo,hi,npt);
1149 }
1150
1152 return curve(start,end,radius,el2,npt,ipt);
1153 }
1154
1155 };
1156
1157
1158
1159 // plot along a line
1160 template<size_t NDIM>
1162 FILE *f=0;
1163 const int npt=traj.npt;
1164
1165 const bool psdot=false;
1166
1167 if(world.rank() == 0) {
1168 f = fopen(filename.c_str(), "w");
1169 if(!f) MADNESS_EXCEPTION("plot_along: failed to open the plot file", 0);
1170
1171 if (psdot) {
1172
1173 fprintf(f,"\\psset{xunit=0.1cm}\n");
1174 fprintf(f,"\\psset{yunit=10cm}\n");
1175 fprintf(f,"\\begin{pspicture}(0,-0.3)(100,1.0)\n");
1176 fprintf(f,"\\pslinewidth=0.05pt\n");
1177 }
1178
1179 // walk along the line
1180 for (int ipt=0; ipt<npt; ipt++) {
1181 Vector<double,NDIM> coord=traj(ipt);
1182 if (psdot) {
1183 long rank=function.evalR(coord);
1185 } else {
1186 fprintf(f,"%4i %12.6f\n",ipt, function(coord));
1187 }
1188 }
1189
1190
1191 if (psdot) fprintf(f,"\\end{pspicture}\n");
1192
1193 fclose(f);
1194 }
1195 world.gop.fence();
1196 ;
1197 }
1198
1199
1200 // plot along a line
1201 template<size_t NDIM>
1202 void plot_along(World& world, trajectory<NDIM> traj, double (*ff)(const Vector<double,NDIM>&), std::string filename) {
1203 FILE *f=0;
1204 const int npt=traj.npt;
1205
1206 const bool psdot=false;
1207
1208 if(world.rank() == 0) {
1209 f = fopen(filename.c_str(), "w");
1210 if(!f) MADNESS_EXCEPTION("plotvtk: failed to open the plot file", 0);
1211
1212 if (psdot) {
1213 fprintf(f,"\\psset{xunit=0.05cm}\n");
1214 fprintf(f,"\\psset{yunit=100cm}\n");
1215 fprintf(f,"\\begin{pspicture}(0,0.25)(100,0.3)\n");
1216 fprintf(f,"\\pslinewidth=0.005pt\n");
1217 }
1218
1219
1220 // walk along the line
1221 for (int ipt=0; ipt<npt; ipt++) {
1222
1223 Vector<double,NDIM> coord=traj(ipt);
1224// fprintf(f,"%12.6f %12.6f %12.6f %12.6f %12.6f %12.6f %12.6f\n",coord[0],coord[1],coord[2],coord[3],coord[4],coord[5], ff(coord));
1225 if (psdot) {
1226 // no hue code here
1227 // long rank=ff.evalR(coord);
1229 } else {
1230 fprintf(f,"%4i %12.6f\n",ipt, ff(coord));
1231 }
1232
1233
1234 }
1235 if (psdot) fprintf(f,"\\end{pspicture}\n");
1236
1237 fclose(f);
1238 }
1239 world.gop.fence();
1240 ;
1241 }
1242
1243
1244
1245}
1246
1247/* @} */
1248#endif // MADNESS_MRA_FUNCPLOT_H__INCLUDED
std::complex< double > double_complex
Definition cfft.h:14
FunctionDefaults holds default paramaters as static class members.
Definition funcdefaults.h:204
static const Tensor< double > & get_cell_width()
Returns the width of each user cell dimension.
Definition funcdefaults.h:468
static const Tensor< double > & get_cell()
Gets the user cell for the simulation.
Definition funcdefaults.h:446
A multiresolution adaptive numerical function.
Definition mra.h:122
Definition indexit.h:180
class for holding the parameters for calculation
Definition QCCalculationParametersBase.h:290
virtual void read_input_and_commandline_options(World &world, const commandlineparser &parser, const std::string tag)
Definition QCCalculationParametersBase.h:325
void set_user_defined_value(const std::string &key, const T &value)
Definition QCCalculationParametersBase.h:533
A tensor is a multidimension array.
Definition tensor.h:317
IsSupported< TensorTypeData< Q >, Tensor< T > & >::type scale(Q x)
Inplace multiplication by scalar of supported type (legacy name)
Definition tensor.h:686
T max(long *ind=0) const
Return the maximum value (and if ind is non-null, its index) in the Tensor.
Definition tensor.h:1703
T min(long *ind=0) const
Return the minimum value (and if ind is non-null, its index) in the Tensor.
Definition tensor.h:1683
A simple, fixed dimension vector.
Definition vector.h:64
size_type size() const
Accessor for the number of elements in the Vector.
Definition vector.h:276
void fence(bool debug=false)
Synchronizes all processes in communicator AND globally ensures no pending AM or tasks.
Definition worldgop.cc:161
void barrier()
Synchronizes all processes in communicator ... does NOT fence pending AM or tasks.
Definition worldgop.h:700
A parallel world class.
Definition world.h:132
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition world.h:318
ProcessID size() const
Returns the number of processes in this World (same as MPI_Comm_size()).
Definition world.h:328
WorldGopInterface & gop
Global operations.
Definition world.h:205
Defines common mathematical and physical constants.
static double lo
Definition dirac-hatom.cc:23
static double function(const coord_3d &r)
Normalized gaussian.
Definition functionio.cc:100
auto T(World &world, response_space &f) -> response_space
Definition global_functions.cc:34
Tensor< typename Tensor< T >::scalar_type > arg(const Tensor< T > &t)
Return a new tensor holding the argument of each element of t (complex types only)
Definition tensor.h:2502
static const double v
Definition hatom_sf_dirac.cc:20
Tensor< double > op(const Tensor< double > &x)
Definition kain.cc:508
#define MADNESS_EXCEPTION(msg, value)
Macro for throwing a MADNESS exception.
Definition madness_exception.h:119
#define MADNESS_ASSERT(condition)
Assert a condition that should be free of side-effects since in release builds this might be a no-op.
Definition madness_exception.h:134
const double pi
Mathematical constant .
Definition constants.h:48
unsigned short htons_x(unsigned short a)
Definition funcplot.h:371
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
static const char * filename
Definition legendre.cc:96
static std::string stringify(T arg)
Definition funcplot.h:1034
Vector< double, 3 > coordT
Definition corepotential.cc:54
response_space scale(response_space a, double b)
void plot_along(World &world, trajectory< NDIM > traj, const Function< double, NDIM > &function, std::string filename)
Definition funcplot.h:1161
std::istream & position_stream_to_word(std::istream &f, const std::string &tag, const char comment='#', bool rewind=true, bool silent=false)
position the input stream to tag, which must be a word (not part of a word)
Definition position_stream.cc:56
void plotvtk_data(const T &function, const char *fieldname, World &world, const char *filename, const Vector< double, NDIM > &plotlo, const Vector< double, NDIM > &plothi, const Vector< long, NDIM > &npt, bool binary=false)
Definition funcplot.h:225
Vector< double, 6 > coord_6d
Definition funcplot.h:1043
void plotdx(const Function< T, NDIM > &f, const char *filename, const Tensor< double > &cell=FunctionDefaults< NDIM >::get_cell(), const std::vector< long > &npt=std::vector< long >(NDIM, 201L), bool binary=true)
Writes an OpenDX format file with a cube/slice of points on a uniform grid.
Definition mraimpl.h:3448
void plot(const std::vector< Function< double, NDIM > > &vf, const std::string &name, const std::vector< std::string > &header)
convenience to get plot_plane and plot_cubefile
Definition funcplot.h:1022
static void plot_line_print_value(FILE *f, double_complex v)
Definition funcplot.h:425
void plot_plane(World &world, const Function< double, NDIM > &function, const std::string name)
Definition funcplot.h:621
void print(const T &t, const Ts &... ts)
Print items to std::cout (items separated by spaces) and terminate with a new line.
Definition print.h:225
NDIM & f
Definition mra.h:2416
NDIM const Function< R, NDIM > & g
Definition mra.h:2416
static void plotpovray(const Function< T, 3 > &function, const char *filename, const Tensor< double > &cell=FunctionDefaults< 3 >::get_cell(), const std::vector< long > &npt=std::vector< long >(3, 201L))
Writes a Povray DF3 format file with a cube of points on a uniform grid.
Definition funcplot.h:382
double imag(double x)
Definition complexfun.h:56
std::enable_if< NDIM==3, void >::type plot_cubefile(World &world, const Function< double, NDIM > &f, std::string filename, std::vector< std::string > molecular_info=std::vector< std::string >(), int npoints=100, double zoom=1.0)
Definition funcplot.h:923
vector< functionT > vecfuncT
Definition corepotential.cc:58
void plotvtk_begin(World &world, const char *filename, const Vector< double, NDIM > &plotlo, const Vector< double, NDIM > &plothi, const Vector< long, NDIM > &npt, bool binary=false)
Definition funcplot.h:146
void plotvtk_end(World &world, const char *filename, bool binary=false)
Definition funcplot.h:349
Vector< double, 3 > coord_3d
Definition funcplot.h:1042
void plot_line(World &world, const char *filename, int npt, const Vector< double, NDIM > &lo, const Vector< double, NDIM > &hi, const opT &op)
Generates ASCII file tabulating f(r) at npoints along line r=lo,...,hi.
Definition funcplot.h:438
double real(double x)
Definition complexfun.h:52
std::string name(const FuncType &type, const int ex=-1)
Definition ccpairfunction.h:28
std::enable_if< NDIM==3, void >::type print_tree_jsonfile(World &world, const Function< double, NDIM > &f, std::string filename)
Definition funcplot.h:998
Function< T, NDIM > copy(const Function< T, NDIM > &f, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &pmap, bool fence=true)
Create a new copy of the function with different distribution and optional fence.
Definition mra.h:2002
static const double b
Definition nonlinschro.cc:119
static const double d
Definition nonlinschro.cc:121
static const double a
Definition nonlinschro.cc:118
static const double c
Definition relops.cc:10
static const double L
Definition rk.cc:46
static const long k
Definition rk.cc:44
Definition test_QCCalculationParametersBase.cc:57
Definition funcplot.h:48
PlotParameters & set_zoom(const double z)
Definition funcplot.h:62
PlotParameters()
Definition funcplot.h:53
std::vector< std::string > plane() const
Definition funcplot.h:93
long npoints() const
Definition funcplot.h:81
double zoom() const
Definition funcplot.h:80
PlotParameters & set_origin(const std::vector< double > origin)
Definition funcplot.h:74
PlotParameters(World &world, const commandlineparser parser=commandlineparser(), const std::string tag="plot")
Definition funcplot.h:49
PlotParameters & set_plane(const std::vector< std::string > plane)
Definition funcplot.h:70
PlotParameters & set_npoints(const long n)
Definition funcplot.h:66
Vector< double, NDIM > origin() const
Definition funcplot.h:84
very simple command line parser
Definition commandlineparser.h:15
Definition funcplot.h:1047
long npt
Definition funcplot.h:1054
static coord_6d circle_6d(const coord_6d &lo, const coord_6d &hi, double radius, coord_3d el2, long npt, long ipt)
Definition funcplot.h:1088
double lo
Definition funcplot.h:1051
static trajectory line_xyz(const int xyz, const long npt)
EZ ctor for a line a direction xyz={0,1,2,..,NDIM-1} thru the origin.
Definition funcplot.h:1143
static trajectory line2(const coordT start, const coordT end, const long npt)
constructor for a line
Definition funcplot.h:1133
trajectory(double radius, coord_3d el2, long npt)
Definition funcplot.h:1120
static Vector< double, NDIM > line_internal(const coordT &lo, const coordT &hi, double radius, coord_3d el2, long npt, long ipt)
Definition funcplot.h:1124
coordT end
Definition funcplot.h:1055
double radius
Definition funcplot.h:1053
coordT start
Definition funcplot.h:1055
static void print_psdot(FILE *f, double x, double y, double color)
Definition funcplot.h:1070
static coord_3d circle2(double lo, double hi, double radius, coord_3d el2, long npt, long ipt)
Definition funcplot.h:1076
trajectory(double radius, long npt)
Definition funcplot.h:1116
trajectory()
Definition funcplot.h:1110
double hi
Definition funcplot.h:1052
Vector< double, NDIM > coordT
Definition funcplot.h:1049
Vector< double, NDIM > operator()(int ipt)
Definition funcplot.h:1151
coord_3d el2
Definition funcplot.h:1056
static double hueCode(const int rank)
some tools for plotting MRA ranks of low order tensors
Definition funcplot.h:1062
coordT(* curve)(const coordT &lo, const coordT &hi, double radius, coord_3d el2, long npt, long ipt)
Definition funcplot.h:1057
InputParameters param
Definition tdse.cc:203
AtomicInt sum
Definition test_atomicint.cc:46
double h(const coord_1d &r)
Definition testgconv.cc:68
static const std::size_t NDIM
Definition testpdiff.cc:42
#define PROFILE_FUNC
Definition worldprofile.h:209