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 
26 
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 
41 namespace Lorene {
42 void 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 
45 using 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 
136 void 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 
165 void 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 
260 void 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 
288 void 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 
380 void 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 
409 void 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 
413 using 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 
507 void 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 }
des_coupe_vtk_y
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
Lorene prototypes.
Definition: app_hor.h:67
Unites
Standard units of space, time and mass.
des_vtk_xyz
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.
des_coupe_vtk_x
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.
des_coupe_vtk_z
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.