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