LORENE
grille_val_interp.C
1 /*
2  * Methods for interpolating with class Grille_val, and its derivative classes.
3  *
4  * See the file grille_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 grille_val_interp_C[] = "$Header: /cvsroot/Lorene/C++/Source/Valencia/grille_val_interp.C,v 1.13 2014/10/13 08:53:48 j_novak Exp $" ;
31 
32 /*
33  * $Id: grille_val_interp.C,v 1.13 2014/10/13 08:53:48 j_novak Exp $
34  * $Log: grille_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 2010/02/04 16:44:35 j_novak
39  * Reformulation of the parabolic interpolation, to have better accuracy
40  *
41  * Revision 1.11 2005/06/23 13:40:08 j_novak
42  * The tests on the number of dimensions have been changed to handle better the
43  * axisymmetric case.
44  *
45  * Revision 1.10 2005/06/22 09:11:17 lm_lin
46  *
47  * Grid wedding: convert from the old C++ object "Cmp" to "Scalar".
48  *
49  * Revision 1.9 2004/05/07 12:32:13 j_novak
50  * New summation from spectral to FD grid. Much faster!
51  *
52  * Revision 1.8 2004/03/25 14:52:33 j_novak
53  * Suppressed some documentation/
54  *
55  * Revision 1.7 2003/12/19 15:05:14 j_novak
56  * Trying to avoid shadowed variables
57  *
58  * Revision 1.6 2003/12/05 14:51:54 j_novak
59  * problem with new SGI compiler
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/09/09 13:00:40 e_gourgoulhon
69  * Modification of declaration of Fortran 77 prototypes for
70  * a better portability (in particular on IBM AIX systems):
71  * All Fortran subroutine names are now written F77_* and are
72  * defined in the new file C++/Include/proto_f77.h.
73  *
74  * Revision 1.2 2001/11/23 16:03:07 j_novak
75  *
76  * minor modifications on the grid check.
77  *
78  * Revision 1.1 2001/11/22 13:41:54 j_novak
79  * Added all source files for manipulating Valencia type objects and making
80  * interpolations to and from Meudon grids.
81  *
82  *
83  * $Header: /cvsroot/Lorene/C++/Source/Valencia/grille_val_interp.C,v 1.13 2014/10/13 08:53:48 j_novak Exp $
84  *
85  */
86 
87 // Fichier includes
88 #include "grille_val.h"
89 #include "proto_f77.h"
90 
91  //------------------
92  // Compatibilite
93  //------------------
94 
95 //Compatibilite entre une grille valencienne cartesienne et une meudonaise
96 namespace Lorene {
97 bool Gval_cart::compatible(const Map* mp, const int lmax, const int lmin)
98  const {
99 
100  //Seulement avec des mappings du genre affine
101  assert( dynamic_cast<const Map_af*>(mp) != 0x0) ;
102 
103  const Mg3d* mgrid = mp->get_mg() ;
104  assert(lmin >= 0 && lmax <= mgrid->get_nzone()) ;
105  int dim_spec = 1 ;
106  for (int i=lmin; i<lmax; i++) {
107  if ((mgrid->get_nt(i) > 1)&&(dim_spec==1)) dim_spec = 2;
108  if (mgrid->get_np(i) > 1) dim_spec = 3;
109  }
110  if (dim_spec != dim.ndim) {
111  cout << "Grille_val::compatibilite: the number of dimensions" << endl ;
112  cout << "of both grids do not coincide!" << endl;
113  abort() ;
114  }
115  if (type_t != mgrid->get_type_t()) {
116  cout << "Grille_val::compatibilite: the symmetries in theta" << endl ;
117  cout << "of both grids do not coincide!" << endl;
118  abort() ;
119  }
120  if (type_p != mgrid->get_type_p()) {
121  cout << "Grille_val::compatibilite: the symmetries in phi" << endl ;
122  cout << "of both grids do not coincide!" << endl;
123  abort() ;
124  }
125 
126  bool dimension = true ;
127  const Coord& rr = mp->r ;
128 
129  double rout = (+rr)(lmax-1, 0, 0, mgrid->get_nr(lmax-1) - 1) ;
130 
131  dimension &= (rout <= *zrmax) ;
132  switch (dim_spec) {
133  case 1:{
134  dimension &= ((+rr)(lmin,0,0,0) >= *zrmin) ;
135  break ;
136  }
137  case 2: {
138  if (mgrid->get_type_t() == SYM)
139  {dimension &= (*zrmin <= 0.) ;}
140  else {
141  dimension &= (*zrmin <= -rout ) ;}
142  dimension &= (*xmin <= 0.) ;
143  dimension &= (*xmax >= rout ) ;
144  break ;
145  }
146  case 3: {
147  if (mgrid->get_type_t() == SYM)
148  {dimension &= (*zrmin <= 0.) ;}
149  else {
150  dimension &= (*zrmin <= -rout) ;}
151  if (mgrid->get_type_p() == SYM) {
152  dimension &= (*ymin <= 0.) ;
153  dimension &= (*xmin <= -rout) ;
154  }
155  else {
156  dimension &= (*xmin <= -rout ) ;
157  dimension &= (*ymin <= -rout ) ;
158  }
159  dimension &= (*xmax >= rout) ;
160  dimension &= (*ymax >= rout) ;
161  break ;
162  }
163  }
164  return dimension ;
165 
166 }
167 //Compatibilite entre une grille valencienne spherique et une meudonaise
168 bool Gval_spher::compatible(const Map* mp, const int lmax, const int lmin)
169  const {
170 
171  //Seulement avec des mappings du genre affine.
172  assert( dynamic_cast<const Map_af*>(mp) != 0x0) ;
173 
174  int dim_spec = 1 ;
175 
176  const Mg3d* mgrid = mp->get_mg() ;
177  for (int i=lmin; i<lmax; i++) {
178  if ((mgrid->get_nt(i) > 1)&&(dim_spec==1)) dim_spec = 2;
179  if (mgrid->get_np(i) > 1) dim_spec = 3;
180  }
181  if (dim_spec > dim.ndim) {
182  cout << "Grille_val::compatibilite: the number of dimensions" << endl ;
183  cout << "of both grids do not coincide!" << endl;
184  cout << "Spectral : " << dim_spec << "D, FD: " << dim.ndim << "D" << endl ;
185  abort() ;
186  }
187  if (type_t != mgrid->get_type_t()) {
188  cout << "Grille_val::compatibilite: the symmetries in theta" << endl ;
189  cout << "of both grids do not coincide!" << endl;
190  abort() ;
191  }
192  if (type_p != mgrid->get_type_p()) {
193  cout << "Grille_val::compatibilite: the symmetries in phi" << endl ;
194  cout << "of both grids do not coincide!" << endl;
195  abort() ;
196  }
197 
198  const Coord& rr = mp->r ;
199 
200  int i_b = mgrid->get_nr(lmax-1) - 1 ;
201 
202  double rmax = (+rr)(lmax-1, 0, 0, i_b) ;
203  double rmin = (+rr)(lmin, 0, 0, 0) ;
204  double valmax = get_zr(dim.dim[0]+nfantome - 1) ;
205  double valmin = get_zr(-nfantome) ;
206 
207  bool dimension = ((rmax <= (valmax)) && (rmin>= (valmin))) ;
208 
209  return dimension ;
210 }
211 
212 // Teste si la grille valencienne cartesienne est contenue dans le mapping
213 // de Meudon (pour le passage Meudon->Valence )
214 bool Gval_cart::contenue_dans(const Map& mp, const int lmax, const int lmin)
215  const {
216  //Seulement avec des mappings du genre affine
217  assert( dynamic_cast<const Map_af*>(&mp) != 0x0) ;
218 
219  const Mg3d* mgrid = mp.get_mg() ;
220  assert(lmin >= 0 && lmax <= mgrid->get_nzone()) ;
221  int dim_spec = 1 ;
222  for (int i=lmin; i<lmax; i++) {
223  if ((mgrid->get_nt(i) > 1)&&(dim_spec==1)) dim_spec = 2;
224  if (mgrid->get_np(i) > 1) dim_spec = 3;
225  }
226  if (dim_spec != dim.ndim) {
227  cout << "Grille_val::contenue_dans: the number of dimensions" << endl ;
228  cout << "of both grids do not coincide!" << endl;
229  abort() ;
230  }
231  if (type_t != mgrid->get_type_t()) {
232  cout << "Grille_val::contenue_dans: the symmetries in theta" << endl ;
233  cout << "of both grids do not coincide!" << endl;
234  abort() ;
235  }
236  if (type_p != mgrid->get_type_p()) {
237  cout << "Grille_val::contenue_dans: the symmetries in phi" << endl ;
238  cout << "of both grids do not coincide!" << endl;
239  abort() ;
240  }
241 
242  bool dimension = true ;
243  const Coord& rr = mp.r ;
244 
245  //For an affine mapping:
246  double radius = (+rr)(lmax-1,0,0,mgrid->get_nr(lmax-1)-1) ;
247  double radius2 = radius*radius ;
248 
249  if (dim_spec ==1) {
250  dimension &= ((+rr)(lmin,0,0,0) <= *zrmin) ;
251  dimension &= (radius >= *zrmax) ;
252  }
253  if (dim_spec ==2) { //## a transformer en switch...
254  dimension &= ((+rr)(lmin,0,0,0)/radius < 1.e-6) ;
255  dimension &= (*xmin >= 0.) ;
256  if (mgrid->get_type_t() == SYM) dimension &= (*zrmin >= 0.) ;
257  double x1 = *xmax ;
258  double z1 = (fabs(*zrmax)>fabs(*zrmin)? *zrmax : *zrmin) ;
259  dimension &= (x1*x1+z1*z1 <= radius2) ;
260  }
261  if (dim_spec == 3) {
262  if (mgrid->get_type_t() == SYM) dimension &= (*zrmin >= 0.) ;
263  if (mgrid->get_type_p() == SYM) dimension &= (*ymin >= 0.) ;
264  double x1 = (fabs(*xmax)>fabs(*xmin)? *xmax : *xmin) ;
265  double y1 = (fabs(*ymax)>fabs(*ymin)? *ymax : *ymin) ;
266  double z1 = (fabs(*zrmax)>fabs(*zrmin)? *zrmax : *zrmin) ;
267  dimension &= (x1*x1+y1*y1+z1*z1 <= radius2) ;
268  }
269  return dimension ;
270 }
271 
272 // Teste si la grille valencienne spherique est contenue dans le mapping
273 // de Meudon (pour le passage Meudon->Valence )
274 bool Gval_spher::contenue_dans(const Map& mp, const int lmax, const int lmin)
275  const {
276 
277  //Seulement avec des mappings du genre affine.
278  assert( dynamic_cast<const Map_af*>(&mp) != 0x0) ;
279 
280  const Mg3d* mgrid = mp.get_mg() ;
281 
282  if (type_t != mgrid->get_type_t()) {
283  cout << "Grille_val::contenue_dans: the symmetries in theta" << endl ;
284  cout << "of both grids do not coincide!" << endl;
285  abort() ;
286  }
287  if (type_p != mgrid->get_type_p()) {
288  cout << "Grille_val::contenue_dans: the symmetries in phi" << endl ;
289  cout << "of both grids do not coincide!" << endl;
290  abort() ;
291  }
292 
293  const Coord& rr = mp.r ;
294 
295  int i_b = mgrid->get_nr(lmax-1) - 1 ;
296 
297  double rmax = (+rr)(lmax-1, 0, 0, i_b) ;
298  double rmin = (+rr)(lmin, 0, 0, 0) ;
299  double valmin = get_zr(0) ;
300  double valmax = get_zr(dim.dim[0] - 1) ;
301 
302  bool dimension = ((rmax >= valmax) && (rmin<= valmin)) ;
303 
304  return dimension ;
305 }
306 
307  //------------------
308  // Interpolation 1D
309  //------------------
310 
311 // Interpolation pour la classe de base
312 Tbl Grille_val::interpol1(const Tbl& rdep, const Tbl& rarr, const Tbl& fdep,
313  int flag, const int type_inter) const {
314  assert(rdep.get_ndim() == 1) ;
315  assert(rarr.get_ndim() == 1) ;
316  assert(rdep.dim == fdep.dim) ;
317 
318  Tbl farr(rarr.dim) ;
319  farr.set_etat_qcq() ;
320 
321  int ndep = rdep.get_dim(0) ;
322  int narr = rarr.get_dim(0) ;
323 
324  switch (type_inter) {
325  case 0: {
326  int ndeg[4] ;
327  ndeg[0] = ndep ;
328  ndeg[1] = narr ;
329  double* err0 = new double[ndep+narr] ;
330  double* err1 = new double[ndep+narr] ;
331  double* den0 = new double[ndep+narr] ;
332  double* den1 = new double[ndep+narr] ;
333  for (int i=0; i<ndep; i++) {
334  err0[i] = rdep(i) ;
335  den0[i] = fdep(i) ;
336  }
337  for (int i=0; i<narr; i++) err1[i] = rarr(i) ;
338  F77_insmts(ndeg, &flag, err0, err1, den0, den1) ;
339  for (int i=0; i<narr; i++) farr.set(i) = den1[i] ;
340  delete[] err0 ;
341  delete[] den0 ;
342  delete[] err1 ;
343  delete[] den1 ;
344  break ;
345  }
346  case 1: {
347  int ip = 0 ;
348  int is = 1 ;
349  assert(ndep > 1);
350  for (int i=0; i<narr; i++) {
351  while(rdep(is) < rarr(i)) is++ ;
352  assert(is<ndep) ;
353  ip = is - 1 ;
354  farr.t[i] = (fdep(is)*(rdep(ip)-rarr(i)) +
355  fdep(ip)*(rarr(i)-rdep(is))) /
356  (rdep(ip)-rdep(is)) ;
357  }
358  break ;
359  }
360 
361  case 2:
362  int i1, i2, i3 ;
363  double xr, x1, x2, x3, y1, y2, y3 ;
364  i2 = 0 ;
365  i3 = 1 ;
366  assert(ndep > 2) ;
367  for (int i=0; i<narr; i++) {
368  xr = rarr(i) ;
369  while(rdep.t[i3] < xr) i3++ ;
370  assert(i3<ndep) ;
371  if (i3 == 1) {
372  i1 = 0 ;
373  i2 = 1 ;
374  i3 = 2 ;
375  }
376  else {
377  i2 = i3 - 1 ;
378  i1 = i2 - 1 ;
379  }
380  x1 = rdep(i1) ;
381  x2 = rdep(i2) ;
382  x3 = rdep(i3) ;
383  y1 = fdep(i1) ;
384  y2 = fdep(i2) ;
385  y3 = fdep(i3) ;
386  double c = y1 ;
387  double b = (y2 - y1) / (x2 - x1) ;
388  double a = ( (y3 - y2)/(x3 - x2) - (y2 - y1)/(x2 - x1) ) / (x3 - x1) ;
389  farr.t[i] = c + b*(xr - x1) + a*(xr - x1)*(xr - x2) ;
390  }
391  break ;
392 
393  case 3:
394  cout << "Spline interpolation not implemented yet!" << endl ;
395  abort() ;
396  break ;
397 
398  default:
399  cout << "Unknown type of interpolation!" << endl ;
400  abort() ;
401  break ;
402  }
403  return farr ;
404 }
405 
406  //------------------
407  // Interpolation 2D
408  //------------------
409 
410 // Interpolation pour les classes derivees
411 Tbl Gval_spher::interpol2(const Tbl& fdep, const Tbl& rarr,
412  const Tbl& tarr, const int type_inter) const
413 {
414  assert(dim.ndim >= 2) ;
415  assert(fdep.get_ndim() == 2) ;
416  assert(rarr.get_ndim() == 1) ;
417  assert(tarr.get_ndim() == 1) ;
418 
419  int ntv = tet->get_dim(0) ;
420  int nrv = zr->get_dim(0) ;
421  int ntm = tarr.get_dim(0) ;
422  int nrm = rarr.get_dim(0) ;
423 
424  Tbl *fdept = new Tbl(nrv) ;
425  fdept->set_etat_qcq() ;
426  Tbl intermediaire(ntv, nrm) ;
427  intermediaire.set_etat_qcq() ;
428 
429  Tbl farr(ntm, nrm) ;
430  farr.set_etat_qcq() ;
431 
432  int job = 1 ;
433  for (int i=0; i<ntv; i++) {
434  for (int j=0; j<nrv; j++) fdept->t[j] = fdep.t[i*nrv+j] ;
435  Tbl fr(interpol1(*zr, rarr, *fdept, job, type_inter)) ;
436  job = 0 ;
437  for (int j=0; j<nrm; j++) intermediaire.t[i*nrm+j] = fr.t[j] ;
438  }
439  delete fdept ;
440 
441  fdept = new Tbl(ntv) ;
442  fdept->set_etat_qcq() ;
443  job = 1 ;
444  for (int i=0; i<nrm; i++) {
445  for (int j=0; j<ntv; j++) fdept->t[j] = intermediaire.t[j*nrm+i] ;
446  Tbl fr(interpol1(*tet, tarr, *fdept, job, type_inter)) ;
447  job = 0 ;
448  for (int j=0; j<ntm; j++) farr.set(j,i) = fr(j) ;
449  }
450  delete fdept ;
451  return farr ;
452 }
453 
454 #ifndef DOXYGEN_SHOULD_SKIP_THIS
455 
456 struct Point {
457  double x ;
458  int l ;
459  int k ;
460 };
461 
462 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
463 
464 int copar(const void* a, const void* b) {
465  double x = (reinterpret_cast<const Point*>(a))->x ;
466  double y = (reinterpret_cast<const Point*>(b))->x ;
467  return x > y ? 1 : -1 ;
468 }
469 
470 Tbl Gval_cart::interpol2(const Tbl& fdep, const Tbl& rarr,
471  const Tbl& tetarr, const int type_inter) const
472 {
473  return interpol2c(*zr, *x, fdep, rarr, tetarr, type_inter) ;
474 }
475 
476 Tbl Gval_cart::interpol2c(const Tbl& zdep, const Tbl& xdep, const Tbl& fdep,
477  const Tbl& rarr, const Tbl& tarr, const int inter_type) const {
478 
479  assert(fdep.get_ndim() == 2) ;
480  assert(zdep.get_ndim() == 1) ;
481  assert(xdep.get_ndim() == 1) ;
482  assert(rarr.get_ndim() == 1) ;
483  assert(tarr.get_ndim() == 1) ;
484 
485  int nz = zdep.get_dim(0) ;
486  int nx = xdep.get_dim(0) ;
487  int nr = rarr.get_dim(0) ;
488  int nt = tarr.get_dim(0) ;
489 
490  Tbl farr(nt, nr) ;
491  farr.set_etat_qcq() ;
492 
493  int narr = nt*nr ;
494  Point* zlk = new Point[narr] ;
495  int inum = 0 ;
496  int ir, it ;
497  for (it=0; it < nt; it++) {
498  for (ir=0; ir < nr; ir++) {
499  zlk[inum].x = rarr(ir)*cos(tarr(it)) ;
500  zlk[inum].l = ir ;
501  zlk[inum].k = it ;
502  inum++ ;
503  }
504  }
505 
506  void* base = reinterpret_cast<void*>(zlk) ;
507  size_t nel = size_t(narr) ;
508  size_t width = sizeof(Point) ;
509  qsort (base, nel, width, copar) ;
510 
511  Tbl effdep(nz) ; effdep.set_etat_qcq() ;
512 
513  double x12 = 1e-6*(zdep(nz-1) - zdep(0)) ;
514  // Attention!! x12 doit etre compatible avec son equivalent dans insmts
515  int ndistz = 0;
516  inum = 0 ;
517  do {
518  inum++ ;
519  if (inum < narr) {
520  if ( (zlk[inum].x - zlk[inum-1].x) > x12 ) {ndistz++ ; }
521  }
522  } while (inum < narr) ;
523  ndistz++ ;
524  Tbl errarr(ndistz) ;
525  errarr.set_etat_qcq() ;
526  Tbl effarr(ndistz) ;
527  ndistz = 0 ;
528  inum = 0 ;
529  do {
530  errarr.set(ndistz) = zlk[inum].x ;
531  inum ++ ;
532  if (inum < narr) {
533  if ( (zlk[inum].x - zlk[inum-1].x) > x12 ) {ndistz++ ; }
534  }
535  } while (inum < narr) ;
536  ndistz++ ;
537 
538  int ijob = 1 ;
539 
540  Tbl tablo(nx, ndistz) ;
541  tablo.set_etat_qcq() ;
542  for (int j=0; j<nx; j++) {
543  for (int i=0; i<nz; i++) effdep.set(i) = fdep(j,i) ;
544  effarr = interpol1(zdep, errarr, effdep, ijob, inter_type) ;
545  ijob = 0 ;
546  for (int i=0; i<ndistz; i++) tablo.set(j,i) = effarr(i) ;
547  }
548 
549  inum = 0 ;
550  int indz = 0 ;
551  Tbl effdep2(nx) ;
552  effdep2.set_etat_qcq() ;
553  while (inum < narr) {
554  Point* xlk = new Point[3*nr] ;
555  int nxline = 0 ;
556  int inum1 ;
557  do {
558  ir = zlk[inum].l ;
559  it = zlk[inum].k ;
560  xlk[nxline].x = rarr(ir)*sin(tarr(it)) ;
561  xlk[nxline].l = ir ;
562  xlk[nxline].k = it ;
563  nxline ++ ; inum ++ ;
564  inum1 = (inum < narr ? inum : 0) ;
565  } while ( ( (zlk[inum1].x - zlk[inum-1].x) < x12 ) && (inum < narr)) ;
566  void* bas2 = reinterpret_cast<void*>(xlk) ;
567  size_t ne2 = size_t(nxline) ;
568  qsort (bas2, ne2, width, copar) ;
569 
570  int inum2 = 0 ;
571  int ndistx = 0 ;
572  do {
573  inum2 ++ ;
574  if (inum2 < nxline) {
575  if ( (xlk[inum2].x - xlk[inum2-1].x) > x12 ) {ndistx++ ; }
576  }
577  } while (inum2 < nxline) ;
578  ndistx++ ;
579 
580  Tbl errarr2(ndistx) ;
581  errarr2.set_etat_qcq() ;
582  Tbl effarr2(ndistx) ;
583  inum2 = 0 ;
584  ndistx = 0 ;
585  do {
586  errarr2.set(ndistx) = xlk[inum2].x ;
587  inum2 ++ ;
588  if (inum2 < nxline) {
589  if ( (xlk[inum2].x - xlk[inum2-1].x) > x12 ) {ndistx++ ; }
590  }
591  } while (inum2 < nxline) ;
592  ndistx++ ;
593 
594  for (int j=0; j<nx; j++) {
595  effdep2.set(j) = tablo(j,indz) ;
596  }
597  indz++ ;
598  ijob = 1 ;
599  effarr2 = interpol1(xdep, errarr2, effdep2, ijob, inter_type) ;
600  int iresu = 0 ;
601  if (ijob == -1) {
602  for (int i=0; i<nxline; i++) {
603  while (fabs(xlk[i].x - xdep(iresu)) > x12 ) {
604  iresu++ ;
605  }
606  ir = xlk[i].l ;
607  it = xlk[i].k ;
608  farr.set(it,ir) = effdep2(iresu) ;
609  }
610  }
611  else {
612  double resu ;
613  for (int i=0; i<nxline; i++) {
614  resu = effarr2(iresu) ;
615  if (i<nxline-1) {
616  if ((xlk[i+1].x-xlk[i].x) > x12) {
617  iresu++ ;
618  }
619  }
620  ir = xlk[i].l ;
621  it = xlk[i].k ;
622  farr.set(it,ir) = resu ;
623  }
624  }
625  delete [] xlk ;
626  }
627 
628  delete [] zlk ;
629  return farr ;
630 }
631 
632 
633  //------------------
634  // Interpolation 3D
635  //------------------
636 
637 // Interpolation pour les classes derivees
638 Tbl Gval_spher::interpol3(const Tbl& fdep, const Tbl& rarr, const Tbl& tarr,
639  const Tbl& parr, const int type_inter) const {
640  assert(dim.ndim == 3) ;
641  assert(fdep.get_ndim() == 3) ;
642  assert(rarr.get_ndim() == 1) ;
643  assert(tarr.get_ndim() == 1) ;
644  assert(parr.get_ndim() == 1) ;
645 
646  int npv = phi->get_dim(0) ;
647  int ntv = tet->get_dim(0) ;
648  int nrv = zr->get_dim(0) ;
649  int npm = parr.get_dim(0) ;
650  int ntm = tarr.get_dim(0) ;
651  int nrm = rarr.get_dim(0) ;
652 
653  Tbl *fdept = new Tbl(ntv, nrv) ;
654  fdept->set_etat_qcq() ;
655  Tbl intermediaire(npv, ntm, nrm) ;
656  intermediaire.set_etat_qcq() ;
657 
658  Tbl farr(npm, ntm, nrm) ;
659  farr.set_etat_qcq() ;
660 
661  for (int i=0; i<npv; i++) {
662  for (int j=0; j<ntv; j++)
663  for (int k=0; k<nrv; k++) fdept->t[j*nrv+k] = fdep.t[(i*ntv+j)*nrv+k] ;
664  Tbl fr(interpol2(*fdept, rarr, tarr, type_inter)) ;
665  for (int j=0; j<ntm; j++)
666  for (int k=0; k<nrm; k++) intermediaire.set(i,j,k) = fr(j,k) ;
667  }
668  delete fdept ;
669 
670  int job = 1 ;
671  fdept = new Tbl(npv) ;
672  fdept->set_etat_qcq() ;
673  for (int i=0; i<ntm; i++) {
674  for (int j=0; j<nrm; j++) {
675  for (int k=0; k<npv; k++) fdept->set(k) = intermediaire(k,i,j) ;
676  Tbl fr(interpol1(*phi, parr, *fdept, job, type_inter)) ;
677  job = 0 ;
678  for (int k=0; k<npm; k++) farr.set(k,i,j) = fr(k) ;
679  }
680  }
681  delete fdept ;
682  return farr ;
683 }
684 
685 Tbl Gval_cart::interpol3(const Tbl& fdep, const Tbl& rarr,
686  const Tbl& tarr, const Tbl& parr, const
687  int inter_type) const {
688 
689  assert(fdep.get_ndim() == 3) ;
690  assert(rarr.get_ndim() == 1) ;
691  assert(tarr.get_ndim() == 1) ;
692  assert(parr.get_ndim() == 1) ;
693 
694  int nz = zr->get_dim(0) ;
695  int nx = x->get_dim(0) ;
696  int ny = y->get_dim(0) ;
697  int nr = rarr.get_dim(0) ;
698  int nt = tarr.get_dim(0) ;
699  int np = parr.get_dim(0) ;
700  Tbl farr(np, nt, nr) ;
701  farr.set_etat_qcq() ;
702 
703  bool coq = (rarr(0)/rarr(nr-1) > 1.e-6) ;
704  Tbl* rarr2(0x0);
705  if (coq) { // If the spectral grid is only made of shells
706  rarr2 = new Tbl(2*nr) ;
707  rarr2->set_etat_qcq() ;
708  double dr = rarr(0)/nr ;
709  for (int i=0; i<nr; i++) rarr2->set(i) = i*dr ;
710  for (int i=nr; i<2*nr; i++) rarr2->set(i) = rarr(i-nr) ;
711  }
712 
713  int nr2 = coq ? 2*nr : nr ;
714 
715  Tbl cylindre(nz, np, nr2) ;
716  cylindre.set_etat_qcq() ;
717  for(int iz=0; iz<nz; iz++) {
718  Tbl carre(ny,nx) ;
719  carre.set_etat_qcq() ;
720  Tbl cercle(np, nr2) ;
721  for (int iy=0; iy<ny; iy++)
722  for (int ix=0; ix<nx; ix++)
723  carre.set(iy,ix) = fdep(iy,ix,iz) ; // This should be optimized...
724  cercle = interpol2c(*x, *y, carre, coq ? *rarr2 : rarr, parr, inter_type) ;
725 
726  for (int ip=0; ip<np; ip++)
727  for (int ir=0; ir<nr2; ir++)
728  cylindre.set(iz,ip,ir) = cercle(ip,ir) ;
729  }
730 
731  for (int ip=0; ip<np; ip++) {
732  Tbl carre(nr2, nz) ;
733  carre.set_etat_qcq() ;
734  Tbl cercle(nt, nr) ;
735  for (int ir=0; ir<nr2; ir++)
736  for (int iz=0; iz<nz; iz++)
737  carre.set(ir,iz) = cylindre(iz,ip,ir) ;
738  cercle = interpol2c(*zr, coq ? *rarr2 : rarr , carre, rarr, tarr,
739  inter_type) ;
740  for (int it=0; it<nt; it++)
741  for (int ir=0; ir<nr; ir++)
742  farr.set(ip,it,ir) = cercle(it,ir) ;
743  }
744 
745  if (coq) delete rarr2 ;
746  return farr ;
747 
748 }
749 
750 
751 }
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::Grille_val::zrmin
double * zrmin
Lower boundary for z (or r ) direction
Definition: grille_val.h:117
Lorene::Grille_val::type_p
int type_p
Type of symmetry in :
Definition: grille_val.h:114
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::Mg3d
Multi-domain grid.
Definition: grilles.h:273
Lorene::Gval_cart::x
Tbl * x
Arrays containing the values of coordinate x on the nodes.
Definition: grille_val.h:343
Lorene::Gval_cart::ymax
double * ymax
Higher boundary for y dimension.
Definition: grille_val.h:339
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::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Gval_cart::interpol3
virtual Tbl interpol3(const Tbl &fdep, const Tbl &rarr, const Tbl &tetarr, const Tbl &phiarr, const int type_inter) const
Performs 3D interpolation.
Definition: grille_val_interp.C:685
Lorene::Map::r
Coord r
r coordinate centered on the grid
Definition: map.h:718
Lorene::Grille_val::zrmax
double * zrmax
Higher boundary for z (or r ) direction
Definition: grille_val.h:120
Lorene::Dim_tbl::dim
int * dim
Array of dimensions (size: ndim).
Definition: dim_tbl.h:102
Lorene::Gval_spher::interpol3
virtual Tbl interpol3(const Tbl &fdep, const Tbl &rarr, const Tbl &tetarr, const Tbl &phiarr, const int type_inter) const
Performs 3D interpolation.
Definition: grille_val_interp.C:638
Lorene::Tbl::get_ndim
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
Definition: tbl.h:400
Lorene::Grille_val::zr
Tbl * zr
Arrays containing the values of coordinate z (or r) on the nodes
Definition: grille_val.h:124
Lorene::Tbl::dim
Dim_tbl dim
Number of dimensions, size,...
Definition: tbl.h:172
Lorene::Gval_spher::phi
Tbl * phi
Arrays containing the values of coordinate on the nodes.
Definition: grille_val.h:591
Lorene::Tbl::t
double * t
The array of double.
Definition: tbl.h:173
Lorene::Gval_spher::compatible
virtual bool compatible(const Map *mp, const int lmax, const int lmin=0) const
Checks if the spectral grid and mapping are compatible with the Grille_val caracteristics for the int...
Definition: grille_val_interp.C:168
Lorene::Gval_cart::xmin
double * xmin
Lower boundary for x dimension.
Definition: grille_val.h:333
Lorene::Grille_val::type_t
int type_t
Type of symmetry in :
Definition: grille_val.h:109
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Lorene::Mg3d::get_type_t
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Definition: grilles.h:485
Lorene::Coord
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
Lorene::Tbl::get_dim
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:403
Lorene::Gval_cart::compatible
virtual bool compatible(const Map *mp, const int lmax, const int lmin=0) const
Checks if the spectral grid and mapping are compatible with the Grille_val caracteristics for the int...
Definition: grille_val_interp.C:97
Lorene::Gval_cart::y
Tbl * y
Arrays containing the values of coordinate y on the nodes.
Definition: grille_val.h:347
Lorene::Gval_cart::xmax
double * xmax
Higher boundary for x dimension.
Definition: grille_val.h:335
Lorene::Gval_cart::interpol2
virtual Tbl interpol2(const Tbl &fdep, const Tbl &rarr, const Tbl &tetarr, const int type_inter) const
Performs 2D interpolation.
Definition: grille_val_interp.C:470
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
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::Map_af
Affine radial mapping.
Definition: map.h:2027
Lorene::Gval_cart::contenue_dans
virtual bool contenue_dans(const Map &mp, const int lmax, const int lmin=0) const
Checks if Gval_cart is contained inside the spectral grid/mapping within the domains [lmin,...
Definition: grille_val_interp.C:214
Lorene::Mg3d::get_type_p
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Definition: grilles.h:495
Lorene::Grille_val::nfantome
int nfantome
The number of hidden cells (same on each side)
Definition: grille_val.h:104
Lorene::Gval_spher::contenue_dans
virtual bool contenue_dans(const Map &mp, const int lmax, const int lmin=0) const
Checks if Gval_spher is contained inside the spectral grid/mapping within the domains [lmin,...
Definition: grille_val_interp.C:274
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
Lorene::Gval_cart::interpol2c
Tbl interpol2c(const Tbl &xdep, const Tbl &zdep, const Tbl &fdep, const Tbl &rarr, const Tbl &tetarr, const int type_inter) const
Same as before, but the coordinates of source points are passed explicitly (xdep, zdep).
Definition: grille_val_interp.C:476
Lorene::Gval_spher::tet
Tbl * tet
Arrays containing the values of coordinate on the nodes.
Definition: grille_val.h:587
Lorene::Dim_tbl::ndim
int ndim
Number of dimensions of the Tbl: can be 1, 2 or 3.
Definition: dim_tbl.h:101
Lorene::Grille_val::dim
Dim_tbl dim
The dimensions of the grid.
Definition: grille_val.h:102
Lorene::Grille_val::get_zr
double get_zr(const int i) const
Read-only of a particular value of the coordinate z (or r ) at the nodes.
Definition: grille_val.h:199
Lorene::Gval_spher::interpol2
virtual Tbl interpol2(const Tbl &fdep, const Tbl &rarr, const Tbl &tetarr, const int type_inter) const
Performs 2D interpolation.
Definition: grille_val_interp.C:411
Lorene::Gval_cart::ymin
double * ymin
Lower boundary for y dimension.
Definition: grille_val.h:337
Lorene::Map
Base class for coordinate mappings.
Definition: map.h:670
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:69