LORENE
des_coupe_vector.C
1/*
2 * Copyright (c) 2005 Eric Gourgoulhon
3 *
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23
24
25/*
26 * $Id: des_coupe_vector.C,v 1.5 2016/12/05 16:18:06 j_novak Exp $
27 * $Log: des_coupe_vector.C,v $
28 * Revision 1.5 2016/12/05 16:18:06 j_novak
29 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
30 *
31 * Revision 1.4 2014/10/13 08:53:22 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.3 2014/10/06 15:16:04 j_novak
35 * Modified #include directives to use c++ syntax.
36 *
37 * Revision 1.2 2008/08/19 06:42:00 j_novak
38 * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
39 * cast-type operations, and constant strings that must be defined as const char*
40 *
41 * Revision 1.1 2005/03/24 22:01:07 e_gourgoulhon
42 * Plot of a vector field represented by a Vector.
43 *
44 *
45 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Graphics/des_coupe_vector.C,v 1.5 2016/12/05 16:18:06 j_novak Exp $
46 *
47 */
48
49
50// Header C
51#include <cmath>
52
53// Header Lorene
54#include "tensor.h"
55#include "graphique.h"
56#include "param.h"
57#include "utilitaires.h"
58#include "unites.h"
59
60namespace Lorene {
61//******************************************************************************
62
63void des_coupe_vect_x(const Vector& vv, double x0, double scale, double sizefl,
64 int nzdes, const char* title, const Scalar* defsurf, double zoom,
65 bool draw_bound, int ny, int nz) {
66
67 const Map& mp = vv.get_mp() ;
68
69 double a1 = mp.val_r(nzdes-1, 1., M_PI/2., 0.) ;
70 double a2 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI/2.) ;
71 double a3 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI) ;
72 double ray = mp.val_r(nzdes-1, 1., 0., 0.) ;
73
74 ray = ( a1 > ray ) ? a1 : ray ;
75 ray = ( a2 > ray ) ? a2 : ray ;
76 ray = ( a3 > ray ) ? a3 : ray ;
77
78 ray *= zoom ;
79
80 double y_min = mp.get_ori_y() - ray ;
81 double y_max = mp.get_ori_y() + ray ;
82 double z_min = mp.get_ori_z() - ray ;
83 double z_max = mp.get_ori_z() + ray ;
84
85 des_coupe_vect_x(vv, x0, scale, sizefl, y_min, y_max, z_min, z_max, title,
86 defsurf, draw_bound, ny, nz) ;
87
88}
89
90
91
92//******************************************************************************
93
94void des_coupe_vect_x(const Vector& vv, double x0, double scale, double
95 sizefl, double y_min, double y_max, double z_min,
96 double z_max, const char* title, const Scalar* defsurf,
97 bool draw_bound, int ny, int nz) {
98
99 using namespace Unites ;
100
101 const Map& mp = vv.get_mp() ;
102
103 if ( vv.get_triad()->identify() != mp.get_bvect_cart().identify() ) {
104 cout <<
105 "des_coupe_vect_x: the vector must be given in Cartesian components !"
106 << endl ;
107 abort() ;
108 }
109
110
111 // Plot of the vector field
112 // ------------------------
113
114 float* vvy = new float[ny*nz] ;
115 float* vvz = new float[ny*nz] ;
116
117 double hy = (y_max - y_min) / double(ny-1) ;
118 double hza = (z_max - z_min) / double(nz-1) ;
119
120 for (int j=0; j<nz; j++) {
121
122 double z = z_min + hza * j ;
123
124 for (int i=0; i<ny; i++) {
125
126 double y = y_min + hy * i ;
127
128 // Computation of (r,theta,phi) :
129 double r, theta, phi ;
130 mp.convert_absolute(x0, y, z, r, theta, phi) ;
131
132 vvy[ny*j+i] = float(vv(2).val_point(r, theta, phi)) ;
133 vvz[ny*j+i] = float(vv(3).val_point(r, theta, phi)) ;
134
135 }
136 }
137
138 float ymin1 = float(y_min / km) ;
139 float ymax1 = float(y_max / km) ;
140 float zmin1 = float(z_min / km) ;
141 float zmax1 = float(z_max / km) ;
142
143 const char* nomy = "y [km]" ;
144 const char* nomz = "z [km]" ;
145
146 if (title == 0x0) {
147 title = "" ;
148 }
149
150 const char* device = 0x0 ;
151 int newgraph = ( (defsurf != 0x0) || draw_bound ) ? 1 : 3 ;
152
153 des_vect(vvy, vvz, ny, nz, ymin1, ymax1, zmin1, zmax1,
154 scale, sizefl, nomy, nomz, title, device, newgraph) ;
155
156
157 delete [] vvy ;
158 delete [] vvz ;
159
160 // Plot of the surface
161 // -------------------
162
163 if (defsurf != 0x0) {
164
165 assert( &(defsurf->get_mp()) == &mp ) ;
166
167 newgraph = draw_bound ? 0 : 2 ;
168
169 des_surface_x(*defsurf, x0, device, newgraph) ;
170
171 } // End of the surface drawing
172
173
174 // Plot of the domains outer boundaries
175 // ------------------------------------
176
177 if (draw_bound) {
178
179 int ndom = mp.get_mg()->get_nzone() ; // total number of domains
180
181 for (int l=0; l<ndom-1; l++) { // loop on the domains (except the
182 // last one)
183
184 newgraph = (l == ndom-2) ? 2 : 0 ;
185
186 des_domaine_x(mp, l, x0, device, newgraph) ;
187 }
188 }
189
190
191}
192
193
194
195//******************************************************************************
196
197void des_coupe_vect_y(const Vector& vv, double y0, double scale, double sizefl,
198 int nzdes, const char* title, const Scalar* defsurf, double zoom,
199 bool draw_bound, int nx, int nz) {
200
201 const Map& mp = vv.get_mp() ;
202
203 double a1 = mp.val_r(nzdes-1, 1., M_PI/2., 0.) ;
204 double a2 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI/2.) ;
205 double a3 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI) ;
206 double ray = mp.val_r(nzdes-1, 1., 0., 0.) ;
207
208 ray = ( a1 > ray ) ? a1 : ray ;
209 ray = ( a2 > ray ) ? a2 : ray ;
210 ray = ( a3 > ray ) ? a3 : ray ;
211
212 ray *= zoom ;
213
214 double x_min = mp.get_ori_x() - ray ;
215 double x_max = mp.get_ori_x() + ray ;
216 double z_min = mp.get_ori_z() - ray ;
217 double z_max = mp.get_ori_z() + ray ;
218
219
220 des_coupe_vect_y(vv, y0, scale, sizefl, x_min, x_max, z_min, z_max, title,
221 defsurf, draw_bound, nx, nz) ;
222
223}
224
225
226
227//******************************************************************************
228
229void des_coupe_vect_y(const Vector& vv, double y0, double scale, double
230 sizefl, double x_min, double x_max, double z_min,
231 double z_max, const char* title, const Scalar* defsurf,
232 bool draw_bound, int nx, int nz) {
233
234 using namespace Unites ;
235
236 const Map& mp = vv.get_mp() ;
237
238 if ( vv.get_triad()->identify() != mp.get_bvect_cart().identify() ) {
239 cout <<
240 "des_coupe_vect_y: the vector must be given in Cartesian components !"
241 << endl ;
242 abort() ;
243 }
244
245
246 // Plot of the vector field
247 // ------------------------
248
249 float* vvx = new float[nx*nz] ;
250 float* vvz = new float[nx*nz] ;
251
252 double hx = (x_max - x_min) / double(nx-1) ;
253 double hza = (z_max - z_min) / double(nz-1) ;
254
255 for (int j=0; j<nz; j++) {
256
257 double z = z_min + hza * j ;
258
259 for (int i=0; i<nx; i++) {
260
261 double x = x_min + hx * i ;
262
263 // Computation of (r,theta,phi) :
264 double r, theta, phi ;
265 mp.convert_absolute(x, y0, z, r, theta, phi) ;
266
267 vvx[nx*j+i] = float(vv(1).val_point(r, theta, phi)) ;
268 vvz[nx*j+i] = float(vv(3).val_point(r, theta, phi)) ;
269
270 }
271 }
272
273 float xmin1 = float(x_min / km) ;
274 float xmax1 = float(x_max / km) ;
275 float zmin1 = float(z_min / km) ;
276 float zmax1 = float(z_max / km) ;
277
278 const char* nomx = "x [km]" ;
279 const char* nomz = "z [km]" ;
280
281 if (title == 0x0) {
282 title = "" ;
283 }
284
285
286 const char* device = 0x0 ;
287 int newgraph = ( (defsurf != 0x0) || draw_bound ) ? 1 : 3 ;
288
289 des_vect(vvx, vvz, nx, nz, xmin1, xmax1, zmin1, zmax1,
290 scale, sizefl, nomx, nomz, title, device, newgraph) ;
291
292
293 delete [] vvx ;
294 delete [] vvz ;
295
296 // Plot of the surface
297 // -------------------
298
299 if (defsurf != 0x0) {
300
301 assert( &(defsurf->get_mp()) == &mp ) ;
302
303 newgraph = draw_bound ? 0 : 2 ;
304
305 des_surface_y(*defsurf, y0, device, newgraph) ;
306
307 } // End of the surface drawing
308
309
310 // Plot of the domains outer boundaries
311 // ------------------------------------
312
313 if (draw_bound) {
314
315 int ndom = mp.get_mg()->get_nzone() ; // total number of domains
316
317 for (int l=0; l<ndom-1; l++) { // loop on the domains (except the
318 // last one)
319
320 newgraph = (l == ndom-2) ? 2 : 0 ;
321
322 des_domaine_y(mp, l, y0, device, newgraph) ;
323 }
324 }
325
326
327}
328
329
330//******************************************************************************
331
332void des_coupe_vect_z(const Vector& vv, double z0, double scale, double sizefl,
333 int nzdes, const char* title, const Scalar* defsurf, double zoom,
334 bool draw_bound, int nx, int ny) {
335
336 const Map& mp = vv.get_mp() ;
337
338 double a1 = mp.val_r(nzdes-1, 1., M_PI/2., 0.) ;
339 double a2 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI/2.) ;
340 double a3 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI) ;
341 double ray = mp.val_r(nzdes-1, 1., 0., 0.) ;
342
343 ray = ( a1 > ray ) ? a1 : ray ;
344 ray = ( a2 > ray ) ? a2 : ray ;
345 ray = ( a3 > ray ) ? a3 : ray ;
346
347 ray *= zoom ;
348
349 double x_min = mp.get_ori_x() - ray ;
350 double x_max = mp.get_ori_x() + ray ;
351 double y_min = mp.get_ori_y() - ray ;
352 double y_max = mp.get_ori_y() + ray ;
353
354 des_coupe_vect_z(vv, z0, scale, sizefl, x_min, x_max, y_min, y_max, title,
355 defsurf, draw_bound, nx, ny) ;
356
357}
358
359
360
361//******************************************************************************
362
363void des_coupe_vect_z(const Vector& vv, double z0, double scale, double
364 sizefl, double x_min, double x_max, double y_min,
365 double y_max, const char* title, const Scalar* defsurf,
366 bool draw_bound, int nx, int ny) {
367
368 using namespace Unites ;
369
370 const Map& mp = vv.get_mp() ;
371
372 if ( vv.get_triad()->identify() != mp.get_bvect_cart().identify() ) {
373 cout <<
374 "des_coupe_vect_y: the vector must be given in Cartesian components !"
375 << endl ;
376 abort() ;
377 }
378
379
380 // Plot of the vector field
381 // ------------------------
382
383 float* vvx = new float[nx*ny] ;
384 float* vvy = new float[nx*ny] ;
385
386 double hy = (y_max - y_min) / double(ny-1) ;
387 double hx = (x_max - x_min) / double(nx-1) ;
388
389 for (int j=0; j<ny; j++) {
390
391 double y = y_min + hy * j ;
392
393 for (int i=0; i<nx; i++) {
394
395 double x = x_min + hx * i ;
396
397 // Computation of (r,theta,phi) :
398 double r, theta, phi ;
399 mp.convert_absolute(x, y, z0, r, theta, phi) ;
400
401 vvx[nx*j+i] = float(vv(1).val_point(r, theta, phi)) ;
402 vvy[nx*j+i] = float(vv(2).val_point(r, theta, phi)) ;
403
404 }
405 }
406
407 float ymin1 = float(y_min / km) ;
408 float ymax1 = float(y_max / km) ;
409 float xmin1 = float(x_min / km) ;
410 float xmax1 = float(x_max / km) ;
411
412 const char* nomy = "y [km]" ;
413 const char* nomx = "x [km]" ;
414
415 if (title == 0x0) {
416 title = "" ;
417 }
418
419 const char* device = 0x0 ;
420 int newgraph = ( (defsurf != 0x0) || draw_bound ) ? 1 : 3 ;
421
422 des_vect(vvx, vvy, nx, ny, xmin1, xmax1, ymin1, ymax1,
423 scale, sizefl, nomx, nomy, title, device, newgraph) ;
424
425
426 delete [] vvx ;
427 delete [] vvy ;
428
429 // Plot of the surface
430 // -------------------
431
432 if (defsurf != 0x0) {
433
434 assert( &(defsurf->get_mp()) == &mp ) ;
435
436 newgraph = draw_bound ? 0 : 2 ;
437
438 des_surface_z(*defsurf, z0, device, newgraph) ;
439
440 } // End of the surface drawing
441
442 // Plot of the domains outer boundaries
443 // ------------------------------------
444
445 if (draw_bound) {
446
447 int ndom = mp.get_mg()->get_nzone() ; // total number of domains
448
449 for (int l=0; l<ndom-1; l++) { // loop on the domains (except the
450 // last one)
451
452 newgraph = (l == ndom-2) ? 2 : 0 ;
453
454 des_domaine_z(mp, l, z0, device, newgraph) ;
455 }
456 }
457
458}
459}
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:393
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene
Lorene prototypes.
Definition app_hor.h:67
Lorene::z
Coord z
z coordinate centered on the grid
Definition map.h:740
Lorene::y
Coord y
y coordinate centered on the grid
Definition map.h:739
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:142
Lorene::phi
Coord phi
coordinate centered on the grid
Definition map.h:732
Lorene::x
Coord x
x coordinate centered on the grid
Definition map.h:738
Lorene::r
Coord r
r coordinate centered on the grid
Definition map.h:730