LORENE
des_vtk_scalar.C
1/*
2 * Copyright (c) 2012-2013 Jason Penner
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 * Written with the assistance of Ian Hawke, University of Southampton
24 */
25
26char des_coupe_vtk_scalar_C[] = "$ $" ;
27
28// Header C
29#include <cmath>
30#include <iostream>
31#include <string>
32
33// Header Lorene
34#include "graphique_vtk.h"
35#include "graphique.h"
36#include "scalar.h"
37#include "param.h"
38#include "utilitaires.h"
39#include "unites.h"
40
41namespace Lorene {
42void des_coupe_vtk_x(const Scalar& uu, double x0, double y_min, double y_max,
43 double z_min, double z_max, const char* title, int ny, int nz) {
44
45using namespace Unites ;
46
47 const Map& mp = uu.get_mp() ;
48
49 int nx = 1;
50 double hx = 0;
51 double x_min = x0;
52
53 int ii=0;
54 int i,j;
55 char c;
56
57 // Dump Scalar Slice to VTK Format
58 // -------------------------------
59
60 double hy = (y_max - y_min) / double(ny-1) ;
61 double hza = (z_max - z_min) / double(nz-1) ;
62
63// output file management
64 char title2[256] = {' '};
65 char title3[256] = {' '};
66
67 while(title[ii]) {
68 c = char(tolower(title[ii]));
69 if(c==' ') {
70 c = '_';
71 }
72 title2[ii] = c ;
73 ii++;
74 }
75 strcat(title3,title2);
76 strcat(title2,"_x.vtk");
77 ofstream myfile;
78 myfile.open (title2);
79
80
81 myfile << "# vtk DataFile Version 3.0" << endl;
82 myfile << "Data produced by mf_evolve" << endl;
83 myfile << "ASCII" << endl;
84 myfile << "DATASET RECTILINEAR_GRID" << endl;
85 myfile << "DIMENSIONS" << " " << nx << " " << ny << " " << nz << endl;
86
87 cout.precision(12);
88 cout.setf(ios::scientific);
89
90 myfile << "X_COORDINATES " << nx << " FLOAT" << endl;
91 for ( i=0; i<nx; i++){
92 myfile << scientific << setprecision(12) << " " << x_min + hx * i << " " ;
93 }
94 myfile << endl;
95
96 myfile << "Y_COORDINATES " << ny << " FLOAT" << endl;
97 for ( i=0; i<ny; i++){
98 myfile << " " << y_min + hy * i << " " ;
99 }
100 myfile << endl;
101
102 myfile << "Z_COORDINATES " << nz << " FLOAT" << endl;
103 for ( i=0; i<nz; i++){
104 myfile << " " << z_min + hza * i << " " ;
105 }
106 myfile << endl;
107
108 myfile << "POINT_DATA " << nx*ny*nz << endl;
109
110 myfile << endl;
111
112 myfile << "SCALARS " << title3 << " FLOAT" << endl;
113
114 myfile << "LOOKUP_TABLE default" << endl;
115
116 for ( j=0; j<nz; j++) {
117
118 double z = z_min + hza * j ;
119
120 for ( i=0; i<ny; i++) {
121
122 double y = y_min + hy * i ;
123
124 // Computation of (r,theta,phi) :
125 double r, theta, phi ;
126 mp.convert_absolute(x0, y, z, r, theta, phi) ;
127 myfile << " " << float(uu.val_point(r, theta, phi)) << endl;
128 }
129 }
130
131 myfile.close();
132}
133
134//******************************************************************************
135
136void des_coupe_vtk_x(const Scalar& uu, double x0, int nzdes, const char* title,
137 double zoom, int ny, int nz) {
138
139 const Map& mp = uu.get_mp() ;
140
141 double a1 = mp.val_r(nzdes-1, 1., M_PI/2., 0.) ;
142 double a2 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI/2.) ;
143 double a3 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI) ;
144 double ray = mp.val_r(nzdes-1, 1., 0., 0.) ;
145
146 ray = ( a1 > ray ) ? a1 : ray ;
147 ray = ( a2 > ray ) ? a2 : ray ;
148 ray = ( a3 > ray ) ? a3 : ray ;
149
150 ray *= zoom ;
151
152 double y_min = mp.get_ori_y() - ray ;
153 double y_max = mp.get_ori_y() + ray ;
154 double z_min = mp.get_ori_z() - ray ;
155 double z_max = mp.get_ori_z() + ray ;
156
157 printf("Printing %s to file\n",title);
158 des_coupe_vtk_x(uu, x0, y_min, y_max, z_min, z_max, title, ny, nz) ;
159
160}
161
162
163//******************************************************************************
164
165void des_coupe_vtk_y(const Scalar& uu, double y0, double x_min, double x_max,
166 double z_min, double z_max, const char* title, int nx, int nz) {
167
168 using namespace Unites ;
169
170 const Map& mp = uu.get_mp() ;
171
172 int ii = 0;
173 int ny = 1;
174 int i,j;
175 char c;
176
177 double y_min = y0;
178 double hy = 0;
179
180 // Plot of isocontours
181 // -------------------
182
183 double hx = (x_max - x_min) / double(nx-1) ;
184 double hza = (z_max - z_min) / double(nz-1) ;
185
186// output file management
187 char title2[256] = {' '};
188 char title3[256] = {' '};
189
190 while(title[ii]) {
191 c = char(tolower(title[ii]));
192 if(c==' ') {
193 c = '_';
194 }
195 title2[ii] = c ;
196 ii++;
197 }
198 strcat(title3,title2);
199 strcat(title2,"_y.vtk");
200 ofstream myfile;
201 myfile.open (title2);
202
203
204 myfile << "# vtk DataFile Version 3.0" << endl;
205 myfile << "Data produced by mf_evolve" << endl;
206 myfile << "ASCII" << endl;
207 myfile << "DATASET RECTILINEAR_GRID" << endl;
208 myfile << "DIMENSIONS" << " " << nx << " " << ny << " " << nz << endl;
209
210 cout.precision(12);
211 cout.setf(ios::scientific);
212
213 myfile << "X_COORDINATES " << nx << " FLOAT" << endl;
214 for ( i=0; i<nx; i++){
215 myfile << scientific << setprecision(12) << " " << x_min + hx * i << " " ;
216 }
217 myfile << endl;
218
219 myfile << "Y_COORDINATES " << ny << " FLOAT" << endl;
220 for ( i=0; i<ny; i++){
221 myfile << " " << y_min + hy * i << " " ;
222 }
223 myfile << endl;
224
225 myfile << "Z_COORDINATES " << nz << " FLOAT" << endl;
226 for ( i=0; i<nz; i++){
227 myfile << " " << z_min + hza * i << " " ;
228 }
229 myfile << endl;
230
231 myfile << "POINT_DATA " << nx*ny*nz << endl;
232
233 myfile << endl;
234
235 myfile << "SCALARS " << title3 << " FLOAT" << endl;
236
237 myfile << "LOOKUP_TABLE default" << endl;
238
239
240 for (j=0; j<nz; j++) {
241
242 double z = z_min + hza * j ;
243
244 for (i=0; i<nx; i++) {
245
246 double x = x_min + hx * i ;
247
248 // Computation of (r,theta,phi) :
249 double r, theta, phi ;
250 mp.convert_absolute(x, y0, z, r, theta, phi) ;
251 myfile << " " << float(uu.val_point(r, theta, phi)) << endl;
252 }
253 }
254
255 myfile.close() ;
256}
257
258//******************************************************************************
259
260void des_coupe_vtk_y(const Scalar& uu, double y0, int nzdes, const char* title,
261 double zoom, int nx, int nz) {
262
263 const Map& mp = uu.get_mp() ;
264
265 double a1 = mp.val_r(nzdes-1, 1., M_PI/2., 0.) ;
266 double a2 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI/2.) ;
267 double a3 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI) ;
268 double ray = mp.val_r(nzdes-1, 1., 0., 0.) ;
269
270 ray = ( a1 > ray ) ? a1 : ray ;
271 ray = ( a2 > ray ) ? a2 : ray ;
272 ray = ( a3 > ray ) ? a3 : ray ;
273
274 ray *= zoom ;
275
276 double x_min = mp.get_ori_x() - ray ;
277 double x_max = mp.get_ori_x() + ray ;
278 double z_min = mp.get_ori_z() - ray ;
279 double z_max = mp.get_ori_z() + ray ;
280
281 printf("Printing %s to file\n",title);
282 des_coupe_vtk_y(uu, y0, x_min, x_max, z_min, z_max, title, nx, nz) ;
283
284}
285
286//******************************************************************************
287
288void des_coupe_vtk_z(const Scalar& uu, double z0, double x_min, double x_max,
289 double y_min, double y_max, const char* title, int nx, int ny) {
290
291 using namespace Unites ;
292
293 const Map& mp = uu.get_mp() ;
294
295 int ii = 0;
296 int nz = 1;
297 int i,j;
298
299 char c;
300
301 double hza = 0;
302 double z_min = z0;
303
304 // Plot of isocontours
305 // -------------------
306
307 double hy = (y_max - y_min) / double(ny-1) ;
308 double hx = (x_max - x_min) / double(nx-1) ;
309
310// output file management
311 char title2[256] = {' '};
312 char title3[256] = {' '};
313
314 while(title[ii]) {
315 c = char(tolower(title[ii]));
316 if(c==' ') {
317 c = '_';
318 }
319 title2[ii] = c ;
320 ii++;
321 }
322 strcat(title3,title2);
323 strcat(title2,"_z.vtk");
324 ofstream myfile;
325 myfile.open (title2);
326
327
328 myfile << "# vtk DataFile Version 3.0" << endl;
329 myfile << "Data produced by mf_evolve" << endl;
330 myfile << "ASCII" << endl;
331 myfile << "DATASET RECTILINEAR_GRID" << endl;
332 myfile << "DIMENSIONS" << " " << nx << " " << ny << " " << nz << endl;
333
334 cout.precision(12);
335 cout.setf(ios::scientific);
336
337 myfile << "X_COORDINATES " << nx << " FLOAT" << endl;
338 for ( i=0; i<nx; i++){
339 myfile << scientific << setprecision(12) << " " << x_min + hx * i << " " ;
340 }
341 myfile << endl;
342
343 myfile << "Y_COORDINATES " << ny << " FLOAT" << endl;
344 for ( i=0; i<ny; i++){
345 myfile << " " << y_min + hy * i << " " ;
346 }
347 myfile << endl;
348
349 myfile << "Z_COORDINATES " << nz << " FLOAT" << endl;
350 for ( i=0; i<nz; i++){
351 myfile << " " << z_min + hza * i << " " ;
352 }
353 myfile << endl;
354
355 myfile << "POINT_DATA " << nx*ny*nz << endl;
356
357 myfile << endl;
358
359 myfile << "SCALARS " << title3 << " FLOAT" << endl;
360
361 myfile << "LOOKUP_TABLE default" << endl;
362 for (j=0; j<ny; j++) {
363
364 double y = y_min + hy * j ;
365
366 for (i=0; i<nx; i++) {
367
368 double x = x_min + hx * i ;
369
370 // Computation of (r,theta,phi) :
371 double r, theta, phi ;
372 mp.convert_absolute(x, y, z0, r, theta, phi) ;
373 myfile << " " << float(uu.val_point(r, theta, phi)) << endl;
374 }
375 }
376 myfile.close();
377}
378/***************************************************************************/
379
380void des_coupe_vtk_z(const Scalar& uu, double z0, int nzdes, const char* title,
381 double zoom, int nx, int ny) {
382
383 const Map& mp = uu.get_mp() ;
384
385 double a1 = mp.val_r(nzdes-1, 1., M_PI/2., 0.) ;
386 double a2 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI/2.) ;
387 double a3 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI) ;
388 double ray = mp.val_r(nzdes-1, 1., 0., 0.) ;
389
390 ray = ( a1 > ray ) ? a1 : ray ;
391 ray = ( a2 > ray ) ? a2 : ray ;
392 ray = ( a3 > ray ) ? a3 : ray ;
393
394 ray *= zoom ;
395
396 double x_min = mp.get_ori_x() - ray ;
397 double x_max = mp.get_ori_x() + ray ;
398 double y_min = mp.get_ori_y() - ray ;
399 double y_max = mp.get_ori_y() + ray ;
400
401 printf("Printing %s to file\n",title);
402 des_coupe_vtk_z(uu, z0, x_min, x_max, y_min, y_max, title, nx, ny) ;
403
404}
405
406//******************************************************************************
407
408
409void des_vtk_xyz(const Scalar& uu, double x_min, double x_max, double y_min,
410 double y_max, double z_min, double z_max,
411 const char* title, int nx, int ny, int nz) {
412
413using namespace Unites ;
414
415 const Map& mp = uu.get_mp() ;
416
417 int ii=0;
418 int i,j,k;
419 char c;
420
421 // Dump Scalar Slice to Ascii
422 // --------------------------
423
424 double hx = (x_max - x_min) / double(nx-1) ;
425 double hy = (y_max - y_min) / double(ny-1) ;
426 double hz = (z_max - z_min) / double(nz-1) ;
427
428// output file management
429 char title2[256] = {' '};
430 char title3[256] = {' '};
431
432 while(title[ii]) {
433 c = char(tolower(title[ii]));
434 if(c==' ') {
435 c = '_';
436 }
437 title2[ii] = c ;
438 ii++;
439 }
440 strcat(title3,title2);
441 strcat(title2,"_xyz.vtk");
442 ofstream myfile;
443 myfile.open (title2);
444
445
446 myfile << "# vtk DataFile Version 3.0" << endl;
447 myfile << "Data produced by mf_evolve" << endl;
448 myfile << "ASCII" << endl;
449 myfile << "DATASET RECTILINEAR_GRID" << endl;
450 myfile << "DIMENSIONS" << " " << nx << " " << ny << " " << nz << endl;
451
452 cout.precision(12);
453 cout.setf(ios::scientific);
454
455 myfile << "X_COORDINATES " << nx << " FLOAT" << endl;
456 for ( i=0; i<nx; i++){
457 myfile << scientific << setprecision(12) << " " << x_min + hx * i << " " ;
458 }
459 myfile << endl;
460
461 myfile << "Y_COORDINATES " << ny << " FLOAT" << endl;
462 for ( i=0; i<ny; i++){
463 myfile << " " << y_min + hy * i << " " ;
464 }
465 myfile << endl;
466
467 myfile << "Z_COORDINATES " << nz << " FLOAT" << endl;
468 for ( i=0; i<nz; i++){
469 myfile << " " << z_min + hz * i << " " ;
470 }
471 myfile << endl;
472
473 myfile << "POINT_DATA " << nx*ny*nz << endl;
474
475 myfile << endl;
476
477 myfile << "SCALARS " << title3 << " FLOAT" << endl;
478
479 myfile << "LOOKUP_TABLE default" << endl;
480
481 for ( k=0; k<nx; k++) {
482
483 double x = x_min + hx * k ;
484
485 for ( j=0; j<nz; j++) {
486
487 double z = z_min + hz * j ;
488
489 for ( i=0; i<ny; i++) {
490
491 double y = y_min + hy * i ;
492
493 // Computation of (r,theta,phi) :
494 double r, theta, phi ;
495 mp.convert_absolute(x, y, z, r, theta, phi) ;
496 // calculate scalar field uu at (r,theta,phi)
497 myfile << " " << uu.val_point(r, theta, phi) << endl;
498 }
499 }
500 }
501
502 myfile.close();
503}
504
505//******************************************************************************
506
507void des_vtk_xyz(const Scalar& uu, int nzdes, const char* title,
508 int nx, int ny, int nz) {
509
510 const Map& mp = uu.get_mp() ;
511
512 double a1 = mp.val_r(nzdes-1, 1., M_PI/2., 0.) ;
513 double a2 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI/2.) ;
514 double a3 = mp.val_r(nzdes-1, 1., M_PI/2., M_PI) ;
515 double ray = mp.val_r(nzdes-1, 1., 0., 0.) ;
516
517 ray = ( a1 > ray ) ? a1 : ray ;
518 ray = ( a2 > ray ) ? a2 : ray ;
519 ray = ( a3 > ray ) ? a3 : ray ;
520
521 double x_min = mp.get_ori_x() - ray ;
522 double x_max = mp.get_ori_x() + ray ;
523 double y_min = mp.get_ori_y() - ray ;
524 double y_max = mp.get_ori_y() + ray ;
525 double z_min = mp.get_ori_z() - ray ;
526 double z_max = mp.get_ori_z() + ray ;
527
528 printf("Printing %s to file\n",title);
529 des_vtk_xyz(uu, x_min, x_max, y_min, y_max, z_min, z_max, title,
530 nx, ny, nz) ;
531
532}
533
534//******************************************************************************
535}
void des_coupe_vtk_x(const Scalar &uu, double x0, int nzdes, const char *title=0x0, double zoom=1.2, int ny=100, int nz=100)
Saves the data for a Scalar in the plane X=constant.
void des_coupe_vtk_z(const Scalar &uu, double z0, int nzdes, const char *title=0x0, double zoom=1.2, int nx=100, int ny=100)
Saves the data for a Scalar in the plane Z=constant.
void des_vtk_xyz(const Scalar &uu, int nzdes, const char *title=0x0, int nx=100, int ny=100, int nz=100)
Saves a three dimensional Scalar.
void des_coupe_vtk_y(const Scalar &uu, double y0, int nzdes, const char *title=0x0, double zoom=1.2, int nx=100, int nz=100)
Saves the data for a Scalar in the plane Y=constant.
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.