KADATH
graph_scalar.cpp
1 /*
2  Copyright 2017 Philippe Grandclement
3 
4  This file is part of Kadath.
5 
6  Kadath 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 3 of the License, or
9  (at your option) any later version.
10 
11  Kadath 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 Kadath. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #include "headcpp.hpp"
21 #include "space.hpp"
22 #include "scalar.hpp"
23 #include "tensor_impl.hpp"
24 
25 
26 //Pgplot stuff
27 extern "C" {
28 int cpgbeg(int unit, const char *file, int nxsub, int nysub);
29 void cpgsch(float size);
30 void cpgslw(int lw);
31 void cpgscf(int font);
32 void cpgenv(float xmin, float xmax, float ymin, float ymax, int just, int axis);
33 void cpglab(const char *xlbl, const char *ylbl, const char *toplbl);
34 void cpgcont(const float *a, int idim, int jdim, int i1, int i2, int j1, int j2, const float *c, int nc, const float *tr);
35 void cpgend(void);
36 }
37 
38 namespace Kadath {
39 void des_equipot(float* uutab, int nx, int ny, float xmin, float xmax,
40  float ymin, float ymax, int ncour, const char* nomx, const char* nomy,
41  const char* title, const char* device = 0x0, int newgraph = 3,
42  int nxpage = 1, int nypage = 1) ;
43 
44 void des_equipot(float* uutab, int nx, int ny, float xmin, float xmax,
45  float ymin, float ymax, int ncour, const char* nomx, const char* nomy,
46  const char* title, const char* device, int newgraph, int nxpage,
47  int nypage) {
48 
49  // Search for the extremal values of the field :
50  // -------------------------------------------
51 
52  float uumin = uutab[0] ;
53  float uumax = uutab[0] ;
54  for (int i=1; i<nx*ny; i++) {
55  uumin = (uutab[i] < uumin) ? uutab[i] : uumin ;
56  uumax = (uutab[i] > uumax) ? uutab[i] : uumax ;
57  }
58 
59  cout << " " << title << " : min, max : " << uumin << " " << uumax
60  << endl ;
61 
62  // Values of equipotentials
63  // -------------------------
64 
65  float* isopot = new float [ncour] ;
66  float hh = float(uumax-uumin) / float(ncour) ;
67  for (int i=0; i<ncour; i++) {
68  isopot[i] = uumin + hh * float(i) ;
69  }
70 
71  // Array defining the grid for pgcont_
72  // -----------------------------------
73  float hx = (xmax - xmin)/float(nx-1) ;
74  float hy = (ymax - ymin)/float(ny-1) ;
75 
76  float tr[6] ;
77  tr[0] = xmin - hx ;
78  tr[1] = hx ;
79  tr[2] = 0 ;
80  tr[3] = ymin - hy ;
81  tr[4] = 0 ;
82  tr[5] = hy ;
83 
84  // Graphics display
85  // ----------------
86 
87  if ( (newgraph == 1) || (newgraph == 3) ) {
88 
89  if (device == 0x0) device = "?" ;
90 
91  int ier = cpgbeg(0, device, nxpage, nypage) ;
92  if (ier != 1) {
93  cout << "des_equipot: problem in opening PGPLOT display !" << endl ;
94  }
95 
96  }
97 
98  // Taille des caracteres:
99  float size = float(1.3) ;
100  cpgsch(size) ;
101 
102  // Epaisseur des traits:
103  int lepais = 1 ;
104  cpgslw(lepais) ;
105 
106  // Fonte axes: caracteres romains:
107  cpgscf(2) ;
108 
109  // Cadre de la figure
110  cpgenv(xmin, xmax, ymin, ymax, 1, 0 ) ;
111  cpglab(nomx,nomy,title) ;
112 
113  // On n'effectue le dessin que si la dynamique est suffisante
114 
115  float dynamique = float(fabs(uumax - uumin)) ;
116 
117  if (dynamique > 1.e-14) {
118 
119  cpgcont(uutab, nx, ny, 1, nx, 1, ny, isopot, ncour, tr) ;
120 
121  }
122 
123  // Closing the graphical output
124  // ----------------------------
125 
126  if ( (newgraph == 2) || (newgraph == 3) ) {
127  cpgend() ;
128  }
129 
130 
131  delete [] isopot ;
132 
133 }
134 
135 
136 void des_coupe (const Scalar& uu, const Point& x0,
137  int num_un, double var_un_min, double var_un_max,
138  int num_deux, double var_deux_min, double var_deux_max,
139  const char* title, const char* axis_one, const char* axis_two, int ncour, int n_un, int n_deux) {
140 
141  assert ((num_un>0) && (num_un<=uu.get_ndim())) ;
142  assert ((num_deux>0) && (num_deux<=uu.get_ndim())) ;
143  assert (num_un != num_deux) ;
144 
145  // Plot of isocontours
146  // -------------------
147  float* uutab = new float[n_un*n_deux] ;
148 
149  double h_un = (var_un_max - var_un_min) / double(n_un-1) ;
150  double h_deux = (var_deux_max - var_deux_min) / double(n_deux-1) ;
151 
152  Point points(x0) ;
153 
154  for (int j=0; j<n_deux; j++) {
155 
156  double var_deux = var_deux_min + h_deux * j ;
157 
158  for (int i=0; i<n_un; i++) {
159 
160  double var_un = var_un_min + h_un * i ;
161 
162  points.set(num_un) = var_un ;
163  points.set(num_deux) = var_deux ;
164 
165  uutab[n_deux*j+i] = float(uu.val_point(points)) ;
166  }
167  }
168 
169  const char* nomy = (axis_two==0x0) ? "" : axis_two ;
170  const char* nomx = (axis_one==0x0) ? "" : axis_one ;
171 
172  const char* titi = (title==0x0) ? "" : title ;
173 
174  char* device = 0x0 ;
175  int newgraph = 3 ;
176 
177  des_equipot(uutab, n_un, n_deux, float(var_un_min), float(var_un_max), float(var_deux_min), float(var_deux_max),
178  ncour, nomx, nomy,
179  titi, device, newgraph) ;
180 
181  delete [] uutab ;
182 }
183 
184 
185 
186 void des_coupe_zeronotdef (const Scalar& uu, const Point& x0,
187  int num_un, double var_un_min, double var_un_max,
188  int num_deux, double var_deux_min, double var_deux_max,
189  const char* title, const char* axis_one, const char* axis_two, int ncour, int n_un, int n_deux) {
190 
191  assert ((num_un>0) && (num_un<=uu.get_ndim())) ;
192  assert ((num_deux>0) && (num_deux<=uu.get_ndim())) ;
193  assert (num_un != num_deux) ;
194 
195  // Plot of isocontours
196  // -------------------
197  float* uutab = new float[n_un*n_deux] ;
198 
199  double h_un = (var_un_max - var_un_min) / double(n_un-1) ;
200  double h_deux = (var_deux_max - var_deux_min) / double(n_deux-1) ;
201 
202  Point points(x0) ;
203 
204  for (int j=0; j<n_deux; j++) {
205 
206  double var_deux = var_deux_min + h_deux * j ;
207 
208  for (int i=0; i<n_un; i++) {
209 
210  double var_un = var_un_min + h_un * i ;
211 
212  points.set(num_un) = var_un ;
213  points.set(num_deux) = var_deux ;
214 
215  uutab[n_deux*j+i] = float(uu.val_point_zeronotdef(points)) ;
216  }
217  }
218 
219  const char* nomy = (axis_two==0x0) ? "" : axis_two ;
220  const char* nomx = (axis_one==0x0) ? "" : axis_one ;
221 
222  const char* titi = (title==0x0) ? "" : title ;
223 
224  char* device = 0x0 ;
225  int newgraph = 3 ;
226 
227  des_equipot(uutab, n_un, n_deux, float(var_un_min), float(var_un_max), float(var_deux_min), float(var_deux_max),
228  ncour, nomx, nomy,
229  titi, device, newgraph) ;
230 
231  delete [] uutab ;
232 }
233 void des_sphere (const Scalar& uu, const Point& x0, double rad, const char* title, int ncour, int n_theta, int n_phi) {
234 
235  // Check dim :
236  if (uu.get_space().get_ndim()!=3) {
237  cerr << "des_sphere only defined for 3-dimensional spaces" << endl ;
238  abort() ;
239  }
240 
241  // Plot of isocontours
242  // -------------------
243  float* uutab = new float[n_theta*n_phi] ;
244 
245  double h_phi = 2*M_PI / double(n_phi-1) ;
246  double h_theta = M_PI / double(n_theta-1) ;
247 
248  double theta = 0 ;
249  double phi ;
250  Point MM(3) ;
251 
252  for (int j=0 ; j<n_theta ; j++) {
253  phi = 0 ;
254  for (int i=0 ; i<n_phi ; i++) {
255 
256 
257  MM.set(1) = rad*sin(theta)*cos(phi) + x0(1) ;
258  MM.set(2) = rad*sin(theta)*sin(phi) + x0(2) ;
259  MM.set(3) = rad*cos(theta) + x0(3) ;
260  uutab[j*n_theta+i] = float(uu.val_point(MM)) ;
261  phi += h_phi ;
262  }
263  theta += h_theta ;
264  }
265 
266  const char* nomy = "theta" ;
267  const char* nomx = "phi" ;
268 
269  const char* titi = (title==0x0) ? "" : title ;
270 
271  char* device = 0x0 ;
272  int newgraph = 3 ;
273 
274  des_equipot(uutab, n_phi, n_theta, 0, float(2*M_PI), 0, float(M_PI), ncour, nomx, nomy, titi, device, newgraph) ;
275 
276  delete [] uutab ;
277 }}