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
28char 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
61namespace Lorene {
62void 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 ] ;
116
117 title_quotes = new char[ strlen(title0)+3 ] ;
118 strcpy(title_quotes, "\"") ;
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 ] ;
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 ] ;
231 strcat(headername, ".dxhead") ;
232 }
233
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 ") ;
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
283void 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}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Base class for pure radial mappings.
Definition map.h:1536
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
int get_dzpuis() const
Returns dzpuis.
Definition scalar.h:557
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Tensor field of valence 1.
Definition vector.h:188
const Scalar & operator()(int) const
Readonly access to a component.
Definition vector.C:305
Vector(const Map &map, int tipe, const Base_vect &triad_i)
Standard constructor.
Definition vector.C:156
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
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tensor.h:303
Lorene prototypes.
Definition app_hor.h:64