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