LORENE
scalar_visu.C
1 /*
2  * 3D visualization of a Scalar via OpenDX
3  *
4  * (see file scalar.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2003 Eric Gourgoulhon
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 
30 /*
31  * $Id: scalar_visu.C,v 1.10 2016/12/05 16:18:19 j_novak Exp $
32  * $Log: scalar_visu.C,v $
33  * Revision 1.10 2016/12/05 16:18:19 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.9 2014/10/13 08:53:47 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.8 2014/10/06 15:16:16 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.7 2005/02/18 13:14:19 j_novak
43  * Changing of malloc/free to new/delete + suppression of some unused variables
44  * (trying to avoid compilation warnings).
45  *
46  * Revision 1.6 2004/03/11 12:07:55 e_gourgoulhon
47  * Added method visu_section_anim.
48  *
49  * Revision 1.5 2003/12/19 15:18:17 j_novak
50  * Shadow variables hunt
51  *
52  * Revision 1.4 2003/12/16 06:32:57 e_gourgoulhon
53  * Added method visu_box.
54  *
55  * Revision 1.3 2003/12/15 08:30:40 p_grandclement
56  * Addition of #include <string.h>
57  *
58  * Revision 1.2 2003/12/14 21:49:14 e_gourgoulhon
59  * Added argument start_dx (to launch OpenDX as a subprocess).
60  *
61  * Revision 1.1 2003/12/11 16:20:25 e_gourgoulhon
62  * First version.
63  *
64  *
65  * $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_visu.C,v 1.10 2016/12/05 16:18:19 j_novak Exp $
66  *
67  */
68 
69 // C headers
70 #include <cstdlib>
71 #include <cstring>
72 
73 // Lorene headers
74 #include "tensor.h"
75 
76  //-----------------------------------------//
77  // visu_section : special cases //
78  //-----------------------------------------//
79 
80 namespace Lorene {
81 void Scalar::visu_section(const char section_type, double aa, double umin,
82  double umax, double vmin, double vmax, const char* title0,
83  const char* filename0, bool start_dx, int nu, int nv) const {
84 
85  Tbl plane(3,3) ;
86  plane.set_etat_qcq() ; // Memory allocation for the Tbl
87 
88  switch (section_type) {
89 
90  case 'x' : {
91  plane.set(0,0) = aa ; // Origin in the plane
92  plane.set(0,1) = 0. ; // (absolute Cartesian coordinates)
93  plane.set(0,2) = 0. ; //
94 
95  plane.set(1,0) = 0. ; // u-coordinate unit vector
96  plane.set(1,1) = 1. ; // (absolute Cartesian components)
97  plane.set(1,2) = 0. ;
98 
99  plane.set(2,0) = 0. ; // v-coordinate unit vector
100  plane.set(2,1) = 0. ; // (absolute Cartesian components)
101  plane.set(2,2) = 1. ;
102 
103  visu_section(plane, umin, umax, vmin, vmax, title0, filename0,
104  start_dx, nu, nv) ;
105  break ;
106  }
107 
108  case 'y' : {
109  plane.set(0,0) = 0. ; // Origin in the plane
110  plane.set(0,1) = aa ; // (absolute Cartesian coordinates)
111  plane.set(0,2) = 0. ; //
112 
113  plane.set(1,0) = 1. ; // u-coordinate unit vector
114  plane.set(1,1) = 0. ; // (absolute Cartesian components)
115  plane.set(1,2) = 0. ;
116 
117  plane.set(2,0) = 0. ; // v-coordinate unit vector
118  plane.set(2,1) = 0. ; // (absolute Cartesian components)
119  plane.set(2,2) = 1. ;
120 
121  visu_section(plane, umin, umax, vmin, vmax, title0, filename0,
122  start_dx, nu, nv) ;
123  break ;
124  }
125 
126  case 'z' : {
127  plane.set(0,0) = 0. ; // Origin in the plane
128  plane.set(0,1) = 0. ; // (absolute Cartesian coordinates)
129  plane.set(0,2) = aa ; //
130 
131  plane.set(1,0) = 1. ; // u-coordinate unit vector
132  plane.set(1,1) = 0. ; // (absolute Cartesian components)
133  plane.set(1,2) = 0. ;
134 
135  plane.set(2,0) = 0. ; // v-coordinate unit vector
136  plane.set(2,1) = 1. ; // (absolute Cartesian components)
137  plane.set(2,2) = 0. ;
138 
139  visu_section(plane, umin, umax, vmin, vmax, title0, filename0,
140  start_dx, nu, nv) ;
141  break ;
142  }
143 
144  default : {
145  cerr << "Scalar::visu_section : unknown type of section ! \n" ;
146  cerr << " section_type = " << section_type << endl ;
147  break ;
148  }
149  }
150 
151 }
152 
153 
154  //-----------------------------------------//
155  // visu_section : general case //
156  //-----------------------------------------//
157 
158 
159 void Scalar::visu_section(const Tbl& plane, double umin, double umax,
160  double vmin, double vmax, const char* title0, const char* filename0,
161  bool start_dx, int nu, int nv) const {
162 
163  char* title ;
164  char* title_quotes ;
165  if (title0 == 0x0) {
166  title = new char[2] ;
167  strcpy(title, " ") ;
168 
169  title_quotes = new char[4] ;
170  strcpy(title_quotes, "\" \"") ;
171  }
172  else {
173  title = new char[ strlen(title0)+1 ] ;
174  strcpy(title, title0) ;
175 
176  title_quotes = new char[ strlen(title0)+3 ] ;
177  strcpy(title_quotes, "\"") ;
178  strcat(title_quotes, title0) ;
179  strcat(title_quotes, "\"") ;
180  }
181 
182  // --------------------------------------------------------
183  // Data file for OpenDX
184  // --------------------------------------------------------
185 
186  char* filename ;
187  if (filename0 == 0x0) {
188  filename = new char[30] ;
189  strcpy(filename, "scalar_section.dxdata") ;
190  }
191  else {
192  filename = new char[ strlen(filename0)+8 ] ;
193  strcpy(filename, filename0) ;
194  strcat(filename, ".dxdata") ;
195  }
196 
197  ofstream fdata(filename) ; // output file
198 
199  fdata << title << "\n" ;
200  fdata << "size : " << nu << " x " << nv << "\n" ;
201  fdata << "u_min = " << umin << " u_max = " << umax << "\n" ;
202  fdata << "v_min = " << vmin << " v_max = " << vmax << "\n" ;
203 
204  // Plane characterization
205  // ----------------------
206 
207  double xa0 = plane(0,0) ;
208  double ya0 = plane(0,1) ;
209  double za0 = plane(0,2) ;
210 
211  double eux = plane(1,0) ;
212  double euy = plane(1,1) ;
213  double euz = plane(1,2) ;
214 
215  double evx = plane(2,0) ;
216  double evy = plane(2,1) ;
217  double evz = plane(2,2) ;
218 
219 
220  // The spectral coefficients are required
221  va.coef() ;
222  const Mtbl_cf& cva = *(va.c_cf) ;
223 
224  // What follows assumes that the mapping is radial:
225  assert( dynamic_cast<const Map_radial*>(mp) != 0x0 ) ;
226 
227  fdata.precision(5) ;
228  fdata.setf(ios::scientific) ;
229 
230  // Loop on the points in the section plane
231  // ---------------------------------------
232  double du = (umax - umin) / double(nu-1) ;
233  double dv = (vmax - vmin) / double(nv-1) ;
234 
235  int npoint = 0 ; // number of data points per line in the file
236 
237  for (int j=0; j<nv; j++) {
238 
239  double v = vmin + dv * j ;
240 
241  for (int i=0; i<nu; i++) {
242 
243  double u = umin + du * i ;
244 
245  double xa = xa0 + u * eux + v * evx ;
246  double ya = ya0 + u * euy + v * evy ;
247  double za = za0 + u * euz + v * evz ;
248 
249 
250  // Values of (r,theta,phi) corresponding to (xa,ya,za) :
251  double rr, th, ph ; // polar coordinates of the mapping associated
252  // to *this
253 
254  mp->convert_absolute(xa, ya, za, rr, th, ph) ;
255 
256  // Values of (l,xi,theta',phi') corresponding to (r,theta,phi):
257  double xi ;
258  int l ;
259 
260  mp->val_lx(rr, th, ph, l, xi) ; // radial mapping assumption
261 
262  // Field value at this point:
263 
264  double ff = cva.val_point(l, xi, th, ph) ;
265 
266  fdata.width(14) ;
267  fdata << ff ;
268  npoint++ ;
269 
270  if (npoint == 6) {
271  fdata << "\n" ;
272  npoint = 0 ;
273  }
274 
275  }
276  }
277 
278  if (npoint != 0) fdata << "\n" ;
279 
280  fdata.close() ;
281 
282  // --------------------------------------------------------
283  // Header file for OpenDX
284  // --------------------------------------------------------
285 
286  char* headername ;
287  if (filename0 == 0x0) {
288  headername = new char[30] ;
289  strcpy(headername, "scalar_section.dxhead") ;
290  }
291  else {
292  headername = new char[ strlen(filename0)+9 ] ;
293  strcpy(headername, filename0) ;
294  strcat(headername, ".dxhead") ;
295  }
296 
297  ofstream fheader(headername) ;
298 
299  fheader << "file = " << filename << endl ;
300  fheader << "grid = " << nu << " x " << nv << endl ;
301  fheader << "format = ascii" << endl ;
302  fheader << "interleaving = record" << endl ;
303  fheader << "majority = column" << endl ;
304  fheader << "header = lines 4" << endl ;
305  fheader << "field = " << title_quotes << endl ;
306  fheader << "structure = scalar" << endl ;
307  fheader << "type = float" << endl ;
308  fheader << "dependency = positions" << endl ;
309  fheader << "positions = regular, regular, " << umin << ", " << du
310  << ", " << vmin << ", " << dv << endl ;
311  fheader << endl ;
312  fheader << "end" << endl ;
313 
314  fheader.close() ;
315 
316 
317  if ( start_dx ) { // Launch of OpenDX
318 
319  char* commande = new char[ strlen(headername) + 60 ] ;
320  strcpy(commande, "ln -s ") ;
321  strcat(commande, headername) ;
322  strcat(commande, " visu_section.dxhead") ;
323 
324  system("rm -f visu_section.dxhead") ;
325  system(commande) ; // ln -s headername visu_section.general
326  system("dx -image visu_section.net &") ;
327 
328  delete [] commande ;
329  }
330 
331  // Final cleaning
332  // --------------
333  delete [] title ;
334  delete [] title_quotes ;
335  delete [] filename ;
336  delete [] headername ;
337 
338 }
339 
340 
341  //------------------------------//
342  // visu_box //
343  //------------------------------//
344 
345 void Scalar::visu_box(double xmin, double xmax, double ymin, double ymax,
346  double zmin, double zmax, const char* title0, const char* filename0,
347  bool start_dx, int nx, int ny, int nz) const {
348 
349  const Scalar* scal ;
350  Scalar* scal_tmp = 0x0 ;
351 
352  // Decrease of dzpuis if dzpuis != 0
353  if ( !check_dzpuis(0) ) {
354  scal_tmp = new Scalar(*this) ;
355  scal_tmp->dec_dzpuis(dzpuis) ;
356  scal = scal_tmp ;
357  }
358  else{
359  scal = this ;
360  }
361 
362  char* title ;
363  char* title_quotes ;
364  if (title0 == 0x0) {
365  title = new char[2] ;
366  strcpy(title, " ") ;
367 
368  title_quotes = new char[4] ;
369  strcpy(title_quotes, "\" \"") ;
370  }
371  else {
372  title = new char[ strlen(title0)+1 ] ;
373  strcpy(title, title0) ;
374 
375  title_quotes = new char[ strlen(title0)+3 ] ;
376  strcpy(title_quotes, "\"") ;
377  strcat(title_quotes, title0) ;
378  strcat(title_quotes, "\"") ;
379  }
380 
381  // --------------------------------------------------------
382  // Data file for OpenDX
383  // --------------------------------------------------------
384 
385  char* filename ;
386  if (filename0 == 0x0) {
387  filename = new char[30] ;
388  strcpy(filename, "scalar_box.dxdata") ;
389  }
390  else {
391  filename = new char[ strlen(filename0)+8 ] ;
392  strcpy(filename, filename0) ;
393  strcat(filename, ".dxdata") ;
394  }
395 
396  ofstream fdata(filename) ; // output file
397 
398  fdata << title << "\n" ;
399  fdata << "size : " << nx << " x " << ny << " x " << nz << "\n" ;
400  fdata << "x_min = " << xmin << " x_max = " << xmax << "\n" ;
401  fdata << "y_min = " << ymin << " y_max = " << ymax << "\n" ;
402  fdata << "z_min = " << zmin << " z_max = " << zmax << "\n" ;
403 
404  // The spectral coefficients are required
405  const Valeur& val = scal->va ;
406  val.coef() ;
407  const Mtbl_cf& cva = *(val.c_cf) ;
408 
409  // What follows assumes that the mapping is radial:
410  assert( dynamic_cast<const Map_radial*>(mp) != 0x0 ) ;
411 
412  fdata.precision(5) ;
413  fdata.setf(ios::scientific) ;
414 
415  // Loop on the points of the drawing box
416  // ---------------------------------------
417  double dx = (xmax - xmin) / double(nx-1) ;
418  double dy = (ymax - ymin) / double(ny-1) ;
419  double dz = (zmax - zmin) / double(nz-1) ;
420 
421  int npoint = 0 ; // number of data points per line in the file
422 
423  for (int k=0; k<nz; k++) {
424 
425  double zz = zmin + dz * k ;
426 
427  for (int j=0; j<ny; j++) {
428 
429  double yy = ymin + dy * j ;
430 
431  for (int i=0; i<nx; i++) {
432 
433  double xx = xmin + dx * i ;
434 
435  // Values of (r,theta,phi) corresponding to (xa,ya,za) :
436  double rr, th, ph ; // polar coordinates of the mapping associated
437  // to *this
438 
439  mp->convert_absolute(xx, yy, zz, rr, th, ph) ;
440 
441  // Values of (l,xi,theta',phi') corresponding to (r,theta,phi):
442  double xi ; int l ;
443 
444  mp->val_lx(rr, th, ph, l, xi) ; // radial mapping assumption
445 
446  // Field value at this point:
447 
448  double ff = cva.val_point(l, xi, th, ph) ;
449 
450  fdata.width(14) ;
451  fdata << ff ;
452  npoint++ ;
453 
454  if (npoint == 9) {
455  fdata << "\n" ;
456  npoint = 0 ;
457  }
458 
459  }
460  }
461 
462  }
463 
464  if (npoint != 0) fdata << "\n" ;
465 
466  fdata.close() ;
467 
468  // --------------------------------------------------------
469  // Header file for OpenDX
470  // --------------------------------------------------------
471 
472  char* headername ;
473  if (filename0 == 0x0) {
474  headername = new char[30] ;
475  strcpy(headername, "scalar_box.dxhead") ;
476  }
477  else {
478  headername = new char[ strlen(filename0)+9 ] ;
479  strcpy(headername, filename0) ;
480  strcat(headername, ".dxhead") ;
481  }
482 
483  ofstream fheader(headername) ;
484 
485  fheader << "file = " << filename << endl ;
486  fheader << "grid = " << nx << " x " << ny << " x " << nz << endl ;
487  fheader << "format = ascii" << endl ;
488  fheader << "interleaving = record" << endl ;
489  fheader << "majority = column" << endl ;
490  fheader << "header = lines 5" << endl ;
491  fheader << "field = " << title_quotes << endl ;
492  fheader << "structure = scalar" << endl ;
493  fheader << "type = float" << endl ;
494  fheader << "dependency = positions" << endl ;
495  fheader << "positions = regular, regular, regular, "
496  << xmin << ", " << dx << ", "
497  << ymin << ", " << dy << ", "
498  << zmin << ", " << dz << endl ;
499  fheader << endl ;
500  fheader << "end" << endl ;
501 
502  fheader.close() ;
503 
504 
505  if ( start_dx ) { // Launch of OpenDX
506 
507  char* commande = new char[ strlen(headername) + 60 ] ;
508  strcpy(commande, "ln -s ") ;
509  strcat(commande, headername) ;
510  strcat(commande, " visu_scalar_box.dxhead") ;
511 
512  system("rm -f visu_scalar_box.dxhead") ;
513  system(commande) ; // ln -s headername visu_section.general
514  system("dx -image visu_scalar_box.net &") ;
515 
516  delete [] commande ;
517  }
518 
519 
520  // Final cleaning
521  // --------------
522 
523  if (scal_tmp != 0x0) delete scal_tmp ;
524  delete [] title ;
525  delete [] title_quotes ;
526  delete [] filename ;
527  delete [] headername ;
528 
529 
530 }
531 
532 
533  //-------------------------------------//
534  // visu_section_anim //
535  //-------------------------------------//
536 
537 
538 void Scalar::visu_section_anim(const char section_type, double aa, double umin,
539  double umax, double vmin, double vmax, int jtime, double ,
540  int jgraph, const char* title, const char* filename_root, bool start_dx,
541  int nu, int nv) const {
542 
543  if ( jtime % jgraph != 0 ) return ;
544 
545  // Preparation of the name of output file
546  // --------------------------------------
547  int k = jtime / jgraph ;
548 
549  char* filename ;
550  if (filename_root == 0x0) {
551  filename = new char[40] ;
552  strcpy(filename, "anim") ;
553  }
554  else {
555  filename = new char[ strlen(filename_root)+10 ] ;
556  strcpy(filename, filename_root) ;
557  }
558 
559  char nomk[5] ;
560  sprintf(nomk, "%04d", k) ;
561  strcat(filename, nomk) ;
562 
563  // Call to visu_section to create the output file
564  // ----------------------------------------------
565 
566  visu_section(section_type, aa, umin, umax, vmin, vmax, title, filename,
567  false, nu, nv) ;
568 
569  // Shall we start OpenDX ?
570  // ---------------------
571 
572  if ( start_dx ) { // Launch of OpenDX
573 
574  system("dx -edit anime.net &") ;
575 
576  }
577 
578  // Final cleaning
579  // --------------
580 
581  delete [] filename ;
582 
583 }
584 }
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:151
Scalar(const Map &mpi)
Constructor from mapping.
Definition: scalar.C:210
Lorene prototypes.
Definition: app_hor.h:67
void visu_section_anim(const char section_type, double aa, double umin, double umax, double vmin, double vmax, int jtime, double ttime, int jgraph=1, const char *title=0x0, const char *filename_root=0x0, bool start_dx=false, int nu=200, int nv=200) const
3D visualization via time evolving plane section (animation).
Definition: scalar_visu.C:538
void visu_box(double xmin, double xmax, double ymin, double ymax, double zmin, double zmax, const char *title0=0x0, const char *filename0=0x0, bool start_dx=true, int nx=40, int ny=40, int nz=40) const
3D visualization (volume rendering) via OpenDX.
Definition: scalar_visu.C:345
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
int dzpuis
Power of r by which the quantity represented by this must be divided in the compactified external do...
Definition: scalar.h:409
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
Valeur va
The numerical value of the Scalar.
Definition: scalar.h:411
void visu_section(const char section_type, double aa, double umin, double umax, double vmin, double vmax, const char *title=0x0, const char *filename=0x0, bool start_dx=true, int nu=200, int nv=200) const
3D visualization via a plane section.
Definition: scalar_visu.C:81
double val_point(int l, double x, double theta, double phi) const
Computes the value of the field represented by *this at an arbitrary point, by means of the spectral ...
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:196
Basic array class.
Definition: tbl.h:164
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition: scalar.C:879
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:301
void convert_absolute(double xx, double yy, double zz, double &rr, double &theta, double &pphi) const
Determines the coordinates corresponding to given absolute Cartesian coordinates (X...
Definition: map.C:305
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
virtual void val_lx(double rr, double theta, double pphi, int &l, double &xi) const =0
Computes the domain index l and the value of corresponding to a point given by its physical coordina...