LORENE
vector_visu.C
1 /*
2  * 3D visualization of a Vector via OpenDX
3  *
4  * (see file vector.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 vector_visu_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/vector_visu.C,v 1.5 2014/10/13 08:53:45 j_novak Exp $" ;
29 
30 /*
31  * $Id: vector_visu.C,v 1.5 2014/10/13 08:53:45 j_novak Exp $
32  * $Log: vector_visu.C,v $
33  * Revision 1.5 2014/10/13 08:53:45 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.4 2014/10/06 15:13:21 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.3 2005/02/16 15:31:56 m_forot
40  * Add the visu_streamline function
41  *
42  * Revision 1.2 2003/12/15 08:30:39 p_grandclement
43  * Addition of #include <string.h>
44  *
45  * Revision 1.1 2003/12/14 21:48:26 e_gourgoulhon
46  * First version
47  *
48  *
49  * $Header: /cvsroot/Lorene/C++/Source/Tensor/vector_visu.C,v 1.5 2014/10/13 08:53:45 j_novak Exp $
50  *
51  */
52 
53 // C headers
54 #include <cstdlib>
55 #include <cstring>
56 
57 // Lorene headers
58 #include "tensor.h"
59 
60 
61 namespace Lorene {
62 void Vector::visu_arrows(double xmin, double xmax, double ymin, double ymax,
63  double zmin, double zmax, const char* title0, const char* filename0,
64  bool start_dx, int nx, int ny, int nz) const {
65 
66  const Vector* vect ;
67  Vector* vect_tmp = 0x0 ;
68 
69  // The drawing is possible only in Cartesian coordinates:
70  if (*triad == mp->get_bvect_cart()) {
71  vect = this ;
72  }
73  else {
74  if (*triad == mp->get_bvect_spher()) {
75  vect_tmp = new Vector(*this) ;
76  vect_tmp->change_triad( mp->get_bvect_cart() ) ;
77  vect = vect_tmp ;
78  }
79  else {
80  cerr << "Vector::visu_arrows : unknown triad !" << endl ;
81  abort() ;
82  }
83  }
84 
85  // The drawing is possible only if dzpuis = 0
86  bool dzpnonzero = false ;
87  for (int i=1; i<=3; i++) {
88  dzpnonzero = dzpnonzero || !(operator()(i).check_dzpuis(0)) ;
89  }
90  if (dzpnonzero) {
91  if (vect_tmp == 0x0) {
92  vect_tmp = new Vector(*this) ;
93  }
94  for (int i=1; i<=3; i++) {
95  Scalar& cvect = vect_tmp->set(i) ;
96  int dzpuis = cvect.get_dzpuis() ;
97  if (dzpuis != 0) {
98  cvect.dec_dzpuis(dzpuis) ;
99  }
100  }
101  vect = vect_tmp ;
102  }
103 
104  char* title ;
105  char* title_quotes ;
106  if (title0 == 0x0) {
107  title = new char[2] ;
108  strcpy(title, " ") ;
109 
110  title_quotes = new char[4] ;
111  strcpy(title_quotes, "\" \"") ;
112  }
113  else {
114  title = new char[ strlen(title0)+1 ] ;
115  strcpy(title, title0) ;
116 
117  title_quotes = new char[ strlen(title0)+3 ] ;
118  strcpy(title_quotes, "\"") ;
119  strcat(title_quotes, title0) ;
120  strcat(title_quotes, "\"") ;
121  }
122 
123  // --------------------------------------------------------
124  // Data file for OpenDX
125  // --------------------------------------------------------
126 
127  char* filename ;
128  if (filename0 == 0x0) {
129  filename = new char[30] ;
130  strcpy(filename, "vector3d_arrows.dxdata") ;
131  }
132  else {
133  filename = new char[ strlen(filename0)+8 ] ;
134  strcpy(filename, filename0) ;
135  strcat(filename, ".dxdata") ;
136  }
137 
138  ofstream fdata(filename) ; // output file
139 
140  fdata << title << "\n" ;
141  fdata << "size : " << nx << " x " << ny << " x " << nz << "\n" ;
142  fdata << "x_min = " << xmin << " x_max = " << xmax << "\n" ;
143  fdata << "y_min = " << ymin << " y_max = " << ymax << "\n" ;
144  fdata << "z_min = " << zmin << " z_max = " << zmax << "\n" ;
145 
146  // The spectral coefficients are required
147  const Valeur& vax = (vect->operator()(1)).get_spectral_va() ;
148  const Valeur& vay = (vect->operator()(2)).get_spectral_va() ;
149  const Valeur& vaz = (vect->operator()(3)).get_spectral_va() ;
150  vax.coef() ;
151  vay.coef() ;
152  vaz.coef() ;
153  const Mtbl_cf& cvax = *(vax.c_cf) ;
154  const Mtbl_cf& cvay = *(vay.c_cf) ;
155  const Mtbl_cf& cvaz = *(vaz.c_cf) ;
156 
157  // What follows assumes that the mapping is radial:
158  assert( dynamic_cast<const Map_radial*>(mp) != 0x0 ) ;
159 
160  fdata.precision(5) ;
161  fdata.setf(ios::scientific) ;
162 
163  // Loop on the points of the drawing box
164  // ---------------------------------------
165  double dx = (xmax - xmin) / double(nx-1) ;
166  double dy = (ymax - ymin) / double(ny-1) ;
167  double dz = (zmax - zmin) / double(nz-1) ;
168 
169  int npoint = 0 ; // number of data points per line in the file
170 
171  for (int k=0; k<nz; k++) {
172 
173  double zz = zmin + dz * k ;
174 
175  for (int j=0; j<ny; j++) {
176 
177  double yy = ymin + dy * j ;
178 
179  for (int i=0; i<nx; i++) {
180 
181  double xx = xmin + dx * i ;
182 
183  // Values of (r,theta,phi) corresponding to (xa,ya,za) :
184  double rr, th, ph ; // polar coordinates of the mapping associated
185  // to *this
186 
187  mp->convert_absolute(xx, yy, zz, rr, th, ph) ;
188 
189  // Values of (l,xi,theta',phi') corresponding to (r,theta,phi):
190  double xi ; int l ;
191 
192  mp->val_lx(rr, th, ph, l, xi) ; // radial mapping assumption
193 
194  // Field value at this point:
195 
196  double vx = cvax.val_point(l, xi, th, ph) ;
197  double vy = cvay.val_point(l, xi, th, ph) ;
198  double vz = cvaz.val_point(l, xi, th, ph) ;
199 
200  fdata.width(14) ; fdata << vx ;
201  fdata.width(14) ; fdata << vy ;
202  fdata.width(14) ; fdata << vz ;
203  npoint++ ;
204 
205  if (npoint == 3) {
206  fdata << "\n" ;
207  npoint = 0 ;
208  }
209 
210  }
211  }
212 
213  }
214 
215  if (npoint != 0) fdata << "\n" ;
216 
217  fdata.close() ;
218 
219  // --------------------------------------------------------
220  // Header file for OpenDX
221  // --------------------------------------------------------
222 
223  char* headername ;
224  if (filename0 == 0x0) {
225  headername = new char[30] ;
226  strcpy(headername, "vector3d_arrows.dxhead") ;
227  }
228  else {
229  headername = new char[ strlen(filename0)+9 ] ;
230  strcpy(headername, filename0) ;
231  strcat(headername, ".dxhead") ;
232  }
233 
234  ofstream fheader(headername) ;
235 
236  fheader << "file = " << filename << endl ;
237  fheader << "grid = " << nx << " x " << ny << " x " << nz << endl ;
238  fheader << "format = ascii" << endl ;
239  fheader << "interleaving = record-vector" << endl ;
240  fheader << "majority = column" << endl ;
241  fheader << "header = lines 5" << endl ;
242  fheader << "field = " << title_quotes << endl ;
243  fheader << "structure = 3-vector" << endl ;
244  fheader << "type = float" << endl ;
245  fheader << "dependency = positions" << endl ;
246  fheader << "positions = regular, regular, regular, "
247  << xmin << ", " << dx << ", "
248  << ymin << ", " << dy << ", "
249  << zmin << ", " << dz << endl ;
250  fheader << endl ;
251  fheader << "end" << endl ;
252 
253  fheader.close() ;
254 
255 
256  if ( start_dx ) { // Launch of OpenDX
257 
258  char* commande = new char[ strlen(headername) + 60 ] ;
259  strcpy(commande, "ln -s ") ;
260  strcat(commande, headername) ;
261  strcat(commande, " visu_vector3d.dxhead") ;
262 
263  system("rm -f visu_vector3d.dxhead") ;
264  system(commande) ; // ln -s headername visu_section.general
265  system("dx -image visu_vector3d.net &") ;
266 
267  delete [] commande ;
268  }
269 
270 
271  // Final cleaning
272  // --------------
273 
274  if (vect_tmp != 0x0) delete vect_tmp ;
275  delete [] title ;
276  delete [] title_quotes ;
277  delete [] filename ;
278  delete [] headername ;
279 
280 
281 }
282 
283 void Vector::visu_streamline(double xmin, double xmax, double ymin, double ymax,
284  double zmin, double zmax, const char* title0, const char* filename0,
285  bool start_dx, int nx, int ny, int nz) const {
286 
287  const Vector* vect ;
288  Vector* vect_tmp = 0x0 ;
289 
290  // The drawing is possible only in Cartesian coordinates:
291  if (*triad == mp->get_bvect_cart()) {
292  vect = this ;
293  }
294  else {
295  if (*triad == mp->get_bvect_spher()) {
296  vect_tmp = new Vector(*this) ;
297  vect_tmp->change_triad( mp->get_bvect_cart() ) ;
298  vect = vect_tmp ;
299  }
300  else {
301  cerr << "Vector::visu_streamline : unknown triad !" << endl ;
302  abort() ;
303  }
304  }
305 
306  // The drawing is possible only if dzpuis = 0
307  bool dzpnonzero = false ;
308  for (int i=1; i<=3; i++) {
309  dzpnonzero = dzpnonzero || !(operator()(i).check_dzpuis(0)) ;
310  }
311  if (dzpnonzero) {
312  if (vect_tmp == 0x0) {
313  vect_tmp = new Vector(*this) ;
314  }
315  for (int i=1; i<=3; i++) {
316  Scalar& cvect = vect_tmp->set(i) ;
317  int dzpuis = cvect.get_dzpuis() ;
318  if (dzpuis != 0) {
319  cvect.dec_dzpuis(dzpuis) ;
320  }
321  }
322  vect = vect_tmp ;
323  }
324 
325  char* title ;
326  char* title_quotes ;
327  if (title0 == 0x0) {
328  title = new char[2] ;
329  strcpy(title, " ") ;
330 
331  title_quotes = new char[4] ;
332  strcpy(title_quotes, "\" \"") ;
333  }
334  else {
335  title = new char[ strlen(title0)+1 ] ;
336  strcpy(title, title0) ;
337 
338  title_quotes = new char[ strlen(title0)+3 ] ;
339  strcpy(title_quotes, "\"") ;
340  strcat(title_quotes, title0) ;
341  strcat(title_quotes, "\"") ;
342  }
343 
344  // --------------------------------------------------------
345  // Data file for OpenDX
346  // --------------------------------------------------------
347 
348  char* filename ;
349  if (filename0 == 0x0) {
350  filename = new char[30] ;
351  strcpy(filename, "vector3d_streamline.dxdata") ;
352  }
353  else {
354  filename = new char[ strlen(filename0)+8 ] ;
355  strcpy(filename, filename0) ;
356  strcat(filename, ".dxdata") ;
357  }
358 
359  ofstream fdata(filename) ; // output file
360 
361  fdata << title << "\n" ;
362  fdata << "size : " << nx << " x " << ny << " x " << nz << "\n" ;
363  fdata << "x_min = " << xmin << " x_max = " << xmax << "\n" ;
364  fdata << "y_min = " << ymin << " y_max = " << ymax << "\n" ;
365  fdata << "z_min = " << zmin << " z_max = " << zmax << "\n" ;
366 
367  // The spectral coefficients are required
368  const Valeur& vax = (vect->operator()(1)).get_spectral_va() ;
369  const Valeur& vay = (vect->operator()(2)).get_spectral_va() ;
370  const Valeur& vaz = (vect->operator()(3)).get_spectral_va() ;
371  vax.coef() ;
372  vay.coef() ;
373  vaz.coef() ;
374  const Mtbl_cf& cvax = *(vax.c_cf) ;
375  const Mtbl_cf& cvay = *(vay.c_cf) ;
376  const Mtbl_cf& cvaz = *(vaz.c_cf) ;
377 
378  // What follows assumes that the mapping is radial:
379  assert( dynamic_cast<const Map_radial*>(mp) != 0x0 ) ;
380 
381  fdata.precision(5) ;
382  fdata.setf(ios::scientific) ;
383 
384  // Loop on the points of the drawing box
385  // ---------------------------------------
386  double dx = (xmax - xmin) / double(nx-1) ;
387  double dy = (ymax - ymin) / double(ny-1) ;
388  double dz = (zmax - zmin) / double(nz-1) ;
389 
390  int npoint = 0 ; // number of data points per line in the file
391 
392  for (int k=0; k<nz; k++) {
393 
394  double zz = zmin + dz * k ;
395 
396  for (int j=0; j<ny; j++) {
397 
398  double yy = ymin + dy * j ;
399 
400  for (int i=0; i<nx; i++) {
401 
402  double xx = xmin + dx * i ;
403 
404  // Values of (r,theta,phi) corresponding to (xa,ya,za) :
405  double rr, th, ph ; // polar coordinates of the mapping associated
406  // to *this
407 
408  mp->convert_absolute(xx, yy, zz, rr, th, ph) ;
409 
410  // Values of (l,xi,theta',phi') corresponding to (r,theta,phi):
411  double xi ; int l ;
412 
413  mp->val_lx(rr, th, ph, l, xi) ; // radial mapping assumption
414 
415  // Field value at this point:
416 
417  double vx = cvax.val_point(l, xi, th, ph) ;
418  double vy = cvay.val_point(l, xi, th, ph) ;
419  double vz = cvaz.val_point(l, xi, th, ph) ;
420 
421  fdata.width(14) ; fdata << vx ;
422  fdata.width(14) ; fdata << vy ;
423  fdata.width(14) ; fdata << vz ;
424  npoint++ ;
425 
426  if (npoint == 3) {
427  fdata << "\n" ;
428  npoint = 0 ;
429  }
430 
431  }
432  }
433 
434  }
435 
436  if (npoint != 0) fdata << "\n" ;
437 
438  fdata.close() ;
439 
440  // --------------------------------------------------------
441  // Header file for OpenDX
442  // --------------------------------------------------------
443 
444  char* headername ;
445  if (filename0 == 0x0) {
446  headername = new char[30] ;
447  strcpy(headername, "vector3d_streamline.dxhead") ;
448  }
449  else {
450  headername = new char[ strlen(filename0)+9 ] ;
451  strcpy(headername, filename0) ;
452  strcat(headername, ".dxhead") ;
453  }
454 
455  ofstream fheader(headername) ;
456 
457  fheader << "file = " << filename << endl ;
458  fheader << "grid = " << nx << " x " << ny << " x " << nz << endl ;
459  fheader << "format = ascii" << endl ;
460  fheader << "interleaving = record-vector" << endl ;
461  fheader << "majority = column" << endl ;
462  fheader << "header = lines 5" << endl ;
463  fheader << "field = " << title_quotes << endl ;
464  fheader << "structure = 3-vector" << endl ;
465  fheader << "type = float" << endl ;
466  fheader << "dependency = positions" << endl ;
467  fheader << "positions = regular, regular, regular, "
468  << xmin << ", " << dx << ", "
469  << ymin << ", " << dy << ", "
470  << zmin << ", " << dz << endl ;
471  fheader << endl ;
472  fheader << "end" << endl ;
473 
474  fheader.close() ;
475 
476 
477  if ( start_dx ) { // Launch of OpenDX
478 
479  char* commande = new char[ strlen(headername) + 60 ] ;
480  strcpy(commande, "ln -s ") ;
481  strcat(commande, headername) ;
482  strcat(commande, " visu_vector3d_SL.general") ;
483 
484  system("rm -f visu_vector3d_SL.general") ;
485  system(commande) ; // ln -s headername visu_section.general
486  system("dx -image visu_vector3d_SL.net &") ;
487 
488  delete [] commande ;
489  }
490 
491 
492  // Final cleaning
493  // --------------
494 
495  if (vect_tmp != 0x0) delete vect_tmp ;
496  delete [] title ;
497  delete [] title_quotes ;
498  delete [] filename ;
499  delete [] headername ;
500 
501 
502 }
503 
504 }
Lorene::Mtbl_cf::val_point
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 ...
Definition: mtbl_cf_val_point.C:115
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
Lorene::Map::val_lx
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...
Lorene::Vector::operator()
const Scalar & operator()(int) const
Readonly access to a component.
Definition: vector.C:305
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Map::get_bvect_spher
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
Lorene::Map_radial
Base class for pure radial mappings.
Definition: map.h:1536
Lorene::Scalar::check_dzpuis
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
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Vector::change_triad
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition: vector_change_triad.C:75
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
Lorene::Scalar::dec_dzpuis
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...
Definition: scalar_r_manip.C:418
Lorene::Vector::Vector
Vector(const Map &map, int tipe, const Base_vect &triad_i)
Standard constructor.
Definition: vector.C:156
Lorene::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:295
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
Lorene::Tensor::triad
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tensor.h:303
Lorene::Map::get_bvect_cart
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition: map.h:791
Lorene::Vector
Tensor field of valence 1.
Definition: vector.h:188
Lorene::Vector::visu_arrows
void visu_arrows(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=8, int ny=8, int nz=8) const
3D visualization via OpenDX.
Definition: vector_visu.C:62
Lorene::Map::convert_absolute
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
Lorene::Scalar::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:557
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:186