LORENE
tbl_val_interp.C
1 /*
2  * Methods for making interpolations with Godunov-type arrays.
3  *
4  * See the file tbl_val.h for documentation
5  *
6  */
7 
8 /*
9  * Copyright (c) 2001 Jerome Novak
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 as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 char TBL_VAL_INTER_C[] = "$Header: /cvsroot/Lorene/C++/Source/Valencia/tbl_val_interp.C,v 1.13 2014/10/13 08:53:48 j_novak Exp $" ;
31 
32 /*
33  * $Id: tbl_val_interp.C,v 1.13 2014/10/13 08:53:48 j_novak Exp $
34  * $Log: tbl_val_interp.C,v $
35  * Revision 1.13 2014/10/13 08:53:48 j_novak
36  * Lorene classes and functions now belong to the namespace Lorene.
37  *
38  * Revision 1.12 2014/10/06 15:13:22 j_novak
39  * Modified #include directives to use c++ syntax.
40  *
41  * Revision 1.11 2013/02/28 16:15:00 j_novak
42  * Minor change.
43  *
44  * Revision 1.10 2007/12/21 10:46:29 j_novak
45  * In "from_spectral..." functions: better treatment of ETATZERO case.
46  *
47  * Revision 1.9 2007/11/02 16:49:12 j_novak
48  * Suppression of intermediate array for spectral summation.
49  *
50  * Revision 1.8 2005/06/23 13:40:08 j_novak
51  * The tests on the number of dimensions have been changed to handle better the
52  * axisymmetric case.
53  *
54  * Revision 1.7 2005/06/22 09:11:17 lm_lin
55  *
56  * Grid wedding: convert from the old C++ object "Cmp" to "Scalar".
57  *
58  * Revision 1.6 2003/12/19 15:05:15 j_novak
59  * Trying to avoid shadowed variables
60  *
61  * Revision 1.5 2003/10/03 16:17:17 j_novak
62  * Corrected some const qualifiers
63  *
64  * Revision 1.4 2002/11/13 11:22:57 j_novak
65  * Version "provisoire" de l'interpolation (sommation depuis la grille
66  * spectrale) aux interfaces de la grille de Valence.
67  *
68  * Revision 1.3 2002/10/16 14:37:15 j_novak
69  * Reorganization of #include instructions of standard C++, in order to
70  * use experimental version 3 of gcc.
71  *
72  * Revision 1.2 2001/11/23 16:03:07 j_novak
73  *
74  * minor modifications on the grid check.
75  *
76  * Revision 1.1 2001/11/22 13:41:54 j_novak
77  * Added all source files for manipulating Valencia type objects and making
78  * interpolations to and from Meudon grids.
79  *
80  *
81  * $Header: /cvsroot/Lorene/C++/Source/Valencia/tbl_val_interp.C,v 1.13 2014/10/13 08:53:48 j_novak Exp $
82  *
83  */
84 
85 // headers C
86 #include <cmath>
87 
88 // headers Lorene
89 #include "headcpp.h"
90 #include "tbl_val.h"
91 
92 
93 namespace Lorene {
94 Scalar Tbl_val::to_spectral(const Map& mp, const int lmax, const int lmin,
95  int type_inter) const {
96 
97  assert(etat != ETATNONDEF) ;
98  assert( gval->compatible(&mp, lmax, lmin) ) ;
99  Scalar resu(mp) ;
100 
101  if (etat == ETATZERO) {
102  resu.annule(lmin, lmax-1) ;
103  return resu ;
104  }
105  else {
106  int nzin = lmax - lmin ;
107  int dim_spec = 1 ;
108  const Mg3d* mgrid = mp.get_mg() ;
109  for (int i=lmin; i<lmax; i++) {
110  if ((mgrid->get_nt(i) > 1)&&(dim_spec==1)) dim_spec = 2;
111  if (mgrid->get_np(i) > 1) dim_spec = 3;
112  }
113  const Coord& rr = mp.r ;
114 
115  int* ntet = new int[nzin] ;
116  int ntetmax = 1 ;
117  int* nphi = new int[nzin] ;
118  int nphimax = 1 ;
119  for (int i=lmin; i<lmax; i++) {
120  int tmp = mgrid->get_np(i) ;
121  nphi[i-lmin] = tmp ;
122  nphimax = (tmp > nphimax ? tmp : nphimax) ;
123  tmp = mgrid->get_nt(i) ;
124  ntet[i-lmin] = tmp ;
125  ntetmax = (tmp > ntetmax ? tmp : ntetmax) ;
126  }
127  if (dim_spec > 1) {
128  for (int i=lmin; i<lmax; i++) {
129  if ((nphimax % nphi[i-lmin]) != 0) {
130  cout << "Tbl_val::to_spectral: The numbers of points in phi" << endl ;
131  cout << "in the different domains of Meudon grid are not" << endl;
132  cout << "well defined; see the documentation." << endl ;
133  abort() ;
134  }
135  assert((ntet[i-lmin]-1) > 0) ;
136  if (((ntetmax-1) % (ntet[i-lmin]-1)) != 0) {
137  cout <<"Tbl_val::to_spectral: The numbers of points in theta"<< endl ;
138  cout << "in the different domains of Meudon grid are not" << endl;
139  cout << "well defined; see the documentation." << endl ;
140  abort() ;
141  }
142  }
143  }
144 
145  resu.allocate_all() ;
146  if (lmin>0) resu.annule(0,lmin-1) ;
147  if (lmax < mgrid->get_nzone()) resu.annule(lmax, mgrid->get_nzone()-1) ;
148 
149  int fant = gval->get_fantome() ;
150  int flag = 1 ;
151  int nrarr = 0 ;
152  for (int l=lmin; l<lmax; l++) nrarr += mgrid->get_nr(l) -1 ;
153  nrarr++ ;
154  switch (dim_spec) {
155 
156  case 1: {
157  int tsize = dim->dim[0] + 2*fant ;
158  Tbl fdep(tsize) ;
159  fdep.set_etat_qcq() ;
160  for (int i=0; i<tsize; i++) fdep.set(i) = t[i] ;
161  Tbl farr(nrarr) ;
162  Tbl rarr(nrarr) ;
163  rarr.set_etat_qcq() ;
164  int inum = 0 ;
165  for (int l=lmin; l<lmax; l++) {
166  for (int i=0; i<mgrid->get_nr(l); i++) {
167  rarr.set(inum) = (+rr)(l,0,0,i) ;
168  inum++ ;
169  }
170  inum--;
171  }
172  farr = gval->interpol1(*gval->zr, rarr, fdep, flag, type_inter) ;
173  inum = 0 ;
174  for (int l=lmin; l<lmax; l++) {
175  for (int i=0; i<mgrid->get_nr(l); i++) {
176  resu.set_grid_point(l,0,0,i) = farr(inum) ;
177  inum++ ;
178  }
179  inum--;
180  }
181  break ;
182  }
183 
184  case 2: {
185  int tsizex = dim->dim[1] + 2*fant ;
186  int tsizez = dim->dim[0] + 2*fant ;
187  Tbl fdep(tsizex, tsizez) ;
188  fdep.set_etat_qcq() ;
189  for (int j=0; j<tsizex; j++) {
190  for (int i=0; i<tsizez; i++) {
191  int l = tsizez*j + i ;
192  fdep.t[l] = t[l] ;
193  }
194  }
195  Tbl farr(ntetmax, nrarr) ;
196  Tbl rarr(nrarr) ;
197  rarr.set_etat_qcq() ;
198  Tbl tetarr(ntetmax) ;
199  tetarr.set_etat_qcq() ;
200  int ltmax = 0 ;
201  int inum = 0 ;
202  for (int l=lmin; l<lmax; l++) {
203  if (ntetmax == ntet[l-lmin]) ltmax = l ;
204  for (int i=0; i<mgrid->get_nr(l); i++) {
205  rarr.set(inum) = (+rr)(l,0,0,i) ;
206  inum++ ;
207  }
208  inum--;
209  }
210  const Coord& tet = mp.tet ;
211  for (int j=0; j<ntetmax; j++)
212  tetarr.set(j) = (+tet)(ltmax,0,j,0) ;
213  farr = gval->interpol2(fdep, rarr, tetarr, type_inter) ;
214  inum = 0 ;
215  for (int l=lmin; l<lmax; l++) {
216  for (int j=0; j<ntet[l-lmin]; j++) {
217  int itet = (ntetmax-1)/(ntet[l-lmin]-1)*j ;
218  for (int i=0; i<mgrid->get_nr(l); i++) {
219  resu.set_grid_point(l,0,j,i) = farr(itet,inum) ;
220  inum++ ;
221  }
222  inum -= mgrid->get_nr(l) ;
223  }
224  inum += mgrid->get_nr(l) - 1;
225  }
226  break ;
227  }
228 
229  case 3: {
230  if (type_inter == 0) {
231  cout << "The use of routine INSMTS is not well suited" << endl ;
232  cout << "for 3D interpolation." << endl ;
233  // abort() ;
234  }
235  int tsizey = dim->dim[2] + 2*fant ;
236  int tsizex = dim->dim[1] + 2*fant ;
237  int tsizez = dim->dim[0] + 2*fant ;
238  Tbl fdep(tsizey, tsizex, tsizez) ;
239  fdep.set_etat_qcq() ;
240  for (int k=0; k<tsizey; k++) {
241  for (int j=0; j<tsizex; j++) {
242  for (int i=0; i<tsizez; i++) {
243  int l = (k*tsizex+j)*tsizez+i ;
244  fdep.t[l] = t[l];
245  }
246  }
247  }
248  Tbl farr(nphimax, ntetmax, nrarr) ;
249  Tbl rarr(nrarr) ;
250  rarr.set_etat_qcq() ;
251  Tbl tetarr(ntetmax) ;
252  tetarr.set_etat_qcq() ;
253  Tbl phiarr(nphimax) ;
254  phiarr.set_etat_qcq() ;
255  int lpmax = 0 ;
256  int ltmax = 0 ;
257  int inum = 0 ;
258  for (int l=lmin; l<lmax; l++) {
259  if (ntetmax == ntet[l-lmin]) ltmax = l ;
260  if (nphimax == nphi[l-lmin]) lpmax = l ;
261  for (int i=0; i<mgrid->get_nr(l); i++) {
262  rarr.set(inum) = (+rr)(l,0,0,i) ;
263  inum++ ;
264  }
265  inum-- ;
266  }
267  const Coord& tet = mp.tet ;
268  const Coord& phi = mp.phi ;
269  for (int k=0; k<nphimax; k++) {
270  phiarr.set(k) = (+phi)(lpmax,k,0,0) ;
271  }
272  for (int j=0; j<ntetmax; j++)
273  tetarr.set(j) = (+tet)(ltmax,0,j,0) ;
274  farr = gval->interpol3(fdep, rarr, tetarr, phiarr, type_inter) ;
275  inum = 0 ;
276  for (int l=lmin; l<lmax; l++) {
277  for (int k=0; k<nphi[l-lmin]; k++) {
278  int iphi = (nphimax-1)/(nphi[l-lmin]-1)*k ;
279  for (int j=0; j<ntet[l-lmin]; j++) {
280  int itet = (ntetmax-1)/(ntet[l-lmin]-1)*j ;
281  for (int i=0; i<mgrid->get_nr(l); i++) {
282  resu.set_grid_point(l,k,j,i) = farr(iphi,itet,inum) ;
283  inum++ ;
284  }
285  inum -= mgrid->get_nr(l) ;
286  }
287  }
288  inum += mgrid->get_nr(l) - 1 ;
289  }
290  break ;
291  }
292 
293  default:
294  cout << "Tbl_val::to_spectral:Strange error..." << endl ;
295  abort() ;
296  break ;
297 
298  }
299 
300  delete [] ntet ;
301  delete [] nphi ;
302  return resu ;
303  }
304 }
305 
306 void Tbl_val::from_spectral(const Scalar& meudon, int lmax, int lmin,
307  bool interfr, bool interft)
308 {
309  assert(meudon.get_etat() != ETATNONDEF) ;
310 #ifndef NDEBUG
311  const Map& mp = meudon.get_mp() ;
312 #endif
313  assert( gval->contenue_dans(mp, lmax, lmin) ) ;
314  if (lmin < 0) {
315  cout << "Tbl_val::from_spectral() : " << endl ;
316  cout << "lmin, lmax : " << lmin << ", " << lmax << endl ;
317  }
318 
319  if (meudon.get_etat() == ETATZERO) {
320  annule_hard() ;
321  return ;
322  }
323  else {
324  assert((meudon.get_etat() == ETATQCQ)||(meudon.get_etat() == ETATUN)) ;
325  set_etat_qcq() ;
326 
327  switch (gval->get_ndim()) {
328 
329  case 1: {
330  gval->somme_spectrale1(meudon, t, get_taille()) ;
331  break ;
332  }
333 
334  case 2: {
335  gval->somme_spectrale2(meudon, t, get_taille()) ;
336  if (interfr) {
337  delete [] tzri ;
338  const Gval_spher* gvs = dynamic_cast<const Gval_spher*>(gval) ; //## A modifier
339  assert (gvs != 0x0) ;
340  tzri = gvs->somme_spectrale2ri(meudon) ;
341  }
342  if (interft) {
343  delete [] txti ;
344  const Gval_spher* gvs = dynamic_cast<const Gval_spher*>(gval) ; //## A modifier
345  assert (gvs != 0x0) ;
346  txti = gvs->somme_spectrale2ti(meudon) ;
347  }
348  break ;
349  }
350 
351  case 3: {
352  gval->somme_spectrale3(meudon, t, get_taille()) ;
353  break ;
354  }
355 
356  default:
357  cout << "Tbl_val::from_spectral:Strange error..." << endl ;
358  abort() ;
359  break ;
360 
361  }
362  return ;
363  }
364 }
365 
366 
367 
368 
369 
370 
371 
372 
373 
374 
375 
376 
377 
378 
379 
380 }
Lorene::Grille_val::interpol1
Tbl interpol1(const Tbl &rdep, const Tbl &rarr, const Tbl &fdep, int flag, const int type_inter) const
Performs 1D interpolation.
Definition: grille_val_interp.C:312
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Lorene::Tbl_val::t
double * t
The array of double at the nodes.
Definition: tbl_val.h:114
Lorene::Grille_val::compatible
virtual bool compatible(const Map *mp, const int lmax, const int lmin=0) const =0
Checks if the spectral grid and mapping are compatible with the Grille_val caracteristics for the int...
Lorene::Tbl_val::etat
int etat
logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition: tbl_val.h:103
Lorene::Tbl_val::from_spectral
void from_spectral(const Scalar &meudon, int lmax, int lmin=0, bool interfr=false, bool interft=false)
Interpolation from a Scalar to a Tbl_val (spectral summation).
Definition: tbl_val_interp.C:306
Lorene::Mg3d
Multi-domain grid.
Definition: grilles.h:273
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Gval_spher::somme_spectrale2ri
double * somme_spectrale2ri(const Scalar &meudon) const
Same as before but at radial grid interfaces.
Definition: gval_from_spectral.C:278
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Grille_val::interpol2
virtual Tbl interpol2(const Tbl &fdep, const Tbl &rarr, const Tbl &tetarr, const int type_inter) const =0
Performs 2D interpolation.
Lorene::Map::r
Coord r
r coordinate centered on the grid
Definition: map.h:718
Lorene::Scalar::set_grid_point
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition: scalar.h:684
Lorene::Grille_val::somme_spectrale3
virtual void somme_spectrale3(const Scalar &meudon, double *t, int taille) const =0
Same as before but for the 3D case.
Lorene::Dim_tbl::dim
int * dim
Array of dimensions (size: ndim).
Definition: dim_tbl.h:102
Lorene::Grille_val::zr
Tbl * zr
Arrays containing the values of coordinate z (or r) on the nodes
Definition: grille_val.h:124
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:861
Lorene::Grille_val::interpol3
virtual Tbl interpol3(const Tbl &fdep, const Tbl &rarr, const Tbl &tetarr, const Tbl &phiarr, const int type_inter) const =0
Performs 3D interpolation.
Lorene::Tbl::t
double * t
The array of double.
Definition: tbl.h:173
Lorene::Tbl_val::tzri
double * tzri
The array at z (or r) interfaces.
Definition: tbl_val.h:116
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Lorene::Coord
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
Lorene::Tbl_val::gval
const Grille_val * gval
The Grille_val (cartesian or spherical) on which the array is defined.
Definition: tbl_val.h:110
Lorene::Tbl_val::to_spectral
Scalar to_spectral(const Map &map, const int lmax, const int lmin=0, int type_inter=2) const
Interpolation from a Tbl_val to a Scalar .
Definition: tbl_val_interp.C:94
Lorene::Grille_val::get_ndim
int get_ndim() const
Returns the number of dimensions.
Definition: grille_val.h:183
Lorene::Tbl_val::annule_hard
void annule_hard()
Sets the Tbl_val to zero in a hard way.
Definition: tbl_val.C:311
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::Grille_val::get_fantome
int get_fantome() const
Returns the number of hidden cells.
Definition: grille_val.h:168
Lorene::Scalar::allocate_all
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition: scalar.C:365
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Map::tet
Coord tet
coordinate centered on the grid
Definition: map.h:719
Lorene::Grille_val::somme_spectrale1
void somme_spectrale1(const Scalar &meudon, double *t, int taille) const
Makes the sommation of the spectral basis functions to know the values of the function described by t...
Definition: gval_from_spectral.C:90
Lorene::Tbl_val::dim
const Dim_tbl * dim
The Dim_tbl giving the dimensions and number of points (without the hidden cells).
Definition: tbl_val.h:108
Lorene::Tbl_val::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl_val.C:294
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Lorene::Grille_val::contenue_dans
virtual bool contenue_dans(const Map &mp, const int lmax, const int lmin=0) const =0
Checks if Grille_val is contained inside the spectral grid/mapping within the domains [lmin,...
Lorene::Map::phi
Coord phi
coordinate centered on the grid
Definition: map.h:720
Lorene::Scalar::annule
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:391
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
Lorene::Tbl_val::get_taille
int get_taille() const
Gives the size of the node array (including the hidden cells)
Definition: tbl_val.h:462
Lorene::Tbl_val::txti
double * txti
The array at x (or ) interfaces.
Definition: tbl_val.h:118
Lorene::Gval_spher::somme_spectrale2ti
double * somme_spectrale2ti(const Scalar &meudon) const
Same as before but at angular grid interfaces.
Definition: gval_from_spectral.C:322
Lorene::Gval_spher
Class for spherical Godunov-type grids.
Definition: grille_val.h:579
Lorene::Map
Base class for coordinate mappings.
Definition: map.h:670
Lorene::Grille_val::somme_spectrale2
virtual void somme_spectrale2(const Scalar &meudon, double *t, int taille) const =0
Same as before but for the 2D case.