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 
46 namespace 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 
150  PROFILE_FUNC;
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 
185  Vector<double, NDIM> space;
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 
243  PROFILE_FUNC;
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 
301  PROFILE_FUNC;
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) {
350  PROFILE_FUNC;
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>
382  static void plotpovray(const Function<T,3>& function,
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) {
440  typedef Vector<double,NDIM> coordT;
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) {
469  typedef Vector<double,NDIM> coordT;
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) {
499  typedef Vector<double,NDIM> coordT;
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) {
532  typedef Vector<double,NDIM> coordT;
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) {
566  typedef Vector<double,NDIM> coordT;
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) {
597  typedef Vector<double,NDIM> coordT;
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>
621  void plot_plane(World& world, const Function<double,NDIM>& function,
622  const std::string name) {
623  typedef std::vector<Function<double,NDIM> > vecfuncT;
624  plot_plane(world,vecfuncT(1,function),name);
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
784 template<size_t NDIM>
785 void 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>
923  plot_cubefile(World& world, const Function<double,NDIM>& f, std::string filename,
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>
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>
1161  void plot_along(World& world, trajectory<NDIM> traj, const Function<double,NDIM>& function, std::string filename) {
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);
1184  trajectory<NDIM>::print_psdot(f,ipt,function(coord),trajectory<NDIM>::hueCode(rank));
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()
Gets the user cell for the simulation.
Definition: funcdefaults.h:446
static const Tensor< double > & get_cell_width()
Returns the width of each user cell dimension.
Definition: funcdefaults.h:468
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
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
IsSupported< TensorTypeData< Q >, Tensor< T > & >::type scale(Q x)
Inplace multiplication by scalar of supported type (legacy name)
Definition: tensor.h:686
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
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 max(a, b)
Definition: lda.h:51
#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
File holds all helper structures necessary for the CC_Operator and CC2 class.
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)
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
void plot_along(World &world, trajectory< NDIM > traj, const Function< double, NDIM > &function, std::string filename)
Definition: funcplot.h:1161
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
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
std::enable_if< NDIM==3, void >::type print_tree_jsonfile(World &world, const Function< double, NDIM > &f, std::string filename)
Definition: funcplot.h:998
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
Function< T, NDIM > sum(World &world, const std::vector< Function< T, NDIM > > &f, bool fence=true)
Returns new function — q = sum_i f[i].
Definition: vmra.h:1421
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::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
std::string type(const PairType &n)
Definition: PNOParameters.h:18
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
static const double b
Definition: nonlinschro.cc:119
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
Vector< double, NDIM > origin() const
Definition: funcplot.h:84
PlotParameters & set_npoints(const long n)
Definition: funcplot.h:66
PlotParameters()
Definition: funcplot.h:53
long npoints() const
Definition: funcplot.h:81
PlotParameters & set_origin(const std::vector< double > origin)
Definition: funcplot.h:74
double zoom() const
Definition: funcplot.h:80
PlotParameters(World &world, const commandlineparser parser=commandlineparser(), const std::string tag="plot")
Definition: funcplot.h:49
PlotParameters & set_zoom(const double z)
Definition: funcplot.h:62
std::vector< std::string > plane() const
Definition: funcplot.h:93
PlotParameters & set_plane(const std::vector< std::string > plane)
Definition: funcplot.h:70
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 Vector< double, NDIM > line_internal(const coordT &lo, const coordT &hi, double radius, coord_3d el2, long npt, long ipt)
Definition: funcplot.h:1124
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
Vector< double, NDIM > operator()(int ipt)
Definition: funcplot.h:1151
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
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
void d()
Definition: test_sig.cc:79
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