30char 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 $" ;
88#include "grille_val.h"
101 assert(
dynamic_cast<const Map_af*
>(mp) != 0x0) ;
104 assert(lmin >= 0 && lmax <= mgrid->get_nzone()) ;
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;
111 cout <<
"Grille_val::compatibilite: the number of dimensions" << endl ;
112 cout <<
"of both grids do not coincide!" << endl;
116 cout <<
"Grille_val::compatibilite: the symmetries in theta" << endl ;
117 cout <<
"of both grids do not coincide!" << endl;
121 cout <<
"Grille_val::compatibilite: the symmetries in phi" << endl ;
122 cout <<
"of both grids do not coincide!" << endl;
126 bool dimension = true ;
129 double rout = (+rr)(lmax-1, 0, 0, mgrid->
get_nr(lmax-1) - 1) ;
131 dimension &= (rout <= *
zrmax) ;
134 dimension &= ((+rr)(lmin,0,0,0) >= *
zrmin) ;
139 {dimension &= (*
zrmin <= 0.) ;}
141 dimension &= (*
zrmin <= -rout ) ;}
142 dimension &= (*
xmin <= 0.) ;
143 dimension &= (*
xmax >= rout ) ;
148 {dimension &= (*
zrmin <= 0.) ;}
150 dimension &= (*
zrmin <= -rout) ;}
152 dimension &= (*
ymin <= 0.) ;
153 dimension &= (*
xmin <= -rout) ;
156 dimension &= (*
xmin <= -rout ) ;
157 dimension &= (*
ymin <= -rout ) ;
159 dimension &= (*
xmax >= rout) ;
160 dimension &= (*
ymax >= rout) ;
172 assert(
dynamic_cast<const Map_af*
>(mp) != 0x0) ;
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;
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 ;
188 cout <<
"Grille_val::compatibilite: the symmetries in theta" << endl ;
189 cout <<
"of both grids do not coincide!" << endl;
193 cout <<
"Grille_val::compatibilite: the symmetries in phi" << endl ;
194 cout <<
"of both grids do not coincide!" << endl;
200 int i_b = mgrid->
get_nr(lmax-1) - 1 ;
202 double rmax = (+rr)(lmax-1, 0, 0, i_b) ;
203 double rmin = (+rr)(lmin, 0, 0, 0) ;
207 bool dimension = ((rmax <= (valmax)) && (rmin>= (valmin))) ;
217 assert(
dynamic_cast<const Map_af*
>(&mp) != 0x0) ;
220 assert(lmin >= 0 && lmax <= mgrid->get_nzone()) ;
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;
227 cout <<
"Grille_val::contenue_dans: the number of dimensions" << endl ;
228 cout <<
"of both grids do not coincide!" << endl;
232 cout <<
"Grille_val::contenue_dans: the symmetries in theta" << endl ;
233 cout <<
"of both grids do not coincide!" << endl;
237 cout <<
"Grille_val::contenue_dans: the symmetries in phi" << endl ;
238 cout <<
"of both grids do not coincide!" << endl;
242 bool dimension = true ;
246 double radius = (+rr)(lmax-1,0,0,mgrid->
get_nr(lmax-1)-1) ;
247 double radius2 = radius*radius ;
250 dimension &= ((+rr)(lmin,0,0,0) <= *
zrmin) ;
251 dimension &= (radius >= *
zrmax) ;
254 dimension &= ((+rr)(lmin,0,0,0)/radius < 1.e-6) ;
255 dimension &= (*
xmin >= 0.) ;
259 dimension &= (x1*x1+z1*z1 <= radius2) ;
267 dimension &= (x1*x1+y1*y1+z1*z1 <= radius2) ;
278 assert(
dynamic_cast<const Map_af*
>(&mp) != 0x0) ;
283 cout <<
"Grille_val::contenue_dans: the symmetries in theta" << endl ;
284 cout <<
"of both grids do not coincide!" << endl;
288 cout <<
"Grille_val::contenue_dans: the symmetries in phi" << endl ;
289 cout <<
"of both grids do not coincide!" << endl;
295 int i_b = mgrid->
get_nr(lmax-1) - 1 ;
297 double rmax = (+rr)(lmax-1, 0, 0, i_b) ;
298 double rmin = (+rr)(lmin, 0, 0, 0) ;
299 double valmin =
get_zr(0) ;
302 bool dimension = ((rmax >= valmax) && (rmin<= valmin)) ;
313 int flag,
const int type_inter)
const {
316 assert(rdep.
dim == fdep.
dim) ;
324 switch (type_inter) {
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++) {
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] ;
350 for (
int i=0; i<narr; i++) {
351 while(rdep(is) < rarr(i)) is++ ;
354 farr.
t[i] = (fdep(is)*(rdep(ip)-rarr(i)) +
355 fdep(ip)*(rarr(i)-rdep(is))) /
356 (rdep(ip)-rdep(is)) ;
363 double xr, x1, x2, x3, y1, y2, y3 ;
367 for (
int i=0; i<narr; i++) {
369 while(rdep.
t[i3] < xr) i3++ ;
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) ;
394 cout <<
"Spline interpolation not implemented yet!" << endl ;
399 cout <<
"Unknown type of interpolation!" << endl ;
412 const Tbl& tarr,
const int type_inter)
const
424 Tbl *fdept =
new Tbl(nrv) ;
426 Tbl intermediaire(ntv, nrm) ;
433 for (
int i=0; i<ntv; i++) {
434 for (
int j=0; j<nrv; j++) fdept->
t[j] = fdep.
t[i*nrv+j] ;
437 for (
int j=0; j<nrm; j++) intermediaire.
t[i*nrm+j] = fr.
t[j] ;
441 fdept =
new Tbl(ntv) ;
444 for (
int i=0; i<nrm; i++) {
445 for (
int j=0; j<ntv; j++) fdept->
t[j] = intermediaire.
t[j*nrm+i] ;
448 for (
int j=0; j<ntm; j++) farr.
set(j,i) = fr(j) ;
454#ifndef DOXYGEN_SHOULD_SKIP_THIS
464int 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 ;
471 const Tbl& tetarr,
const int type_inter)
const
477 const Tbl& rarr,
const Tbl& tarr,
const int inter_type)
const {
494 Point* zlk =
new Point[narr] ;
497 for (it=0; it < nt; it++) {
498 for (ir=0; ir < nr; ir++) {
499 zlk[inum].x = rarr(ir)*
cos(tarr(it)) ;
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) ;
513 double x12 = 1e-6*(zdep(nz-1) - zdep(0)) ;
520 if ( (zlk[inum].
x - zlk[inum-1].
x) > x12 ) {ndistz++ ; }
522 }
while (inum < narr) ;
530 errarr.
set(ndistz) = zlk[inum].x ;
533 if ( (zlk[inum].
x - zlk[inum-1].
x) > x12 ) {ndistz++ ; }
535 }
while (inum < narr) ;
540 Tbl tablo(nx, ndistz) ;
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) ;
546 for (
int i=0; i<ndistz; i++) tablo.
set(j,i) = effarr(i) ;
553 while (inum < narr) {
554 Point* xlk =
new Point[3*nr] ;
560 xlk[nxline].x = rarr(ir)*
sin(tarr(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) ;
574 if (inum2 < nxline) {
575 if ( (xlk[inum2].
x - xlk[inum2-1].
x) > x12 ) {ndistx++ ; }
577 }
while (inum2 < nxline) ;
580 Tbl errarr2(ndistx) ;
582 Tbl effarr2(ndistx) ;
586 errarr2.
set(ndistx) = xlk[inum2].x ;
588 if (inum2 < nxline) {
589 if ( (xlk[inum2].
x - xlk[inum2-1].
x) > x12 ) {ndistx++ ; }
591 }
while (inum2 < nxline) ;
594 for (
int j=0; j<nx; j++) {
595 effdep2.
set(j) = tablo(j,indz) ;
599 effarr2 =
interpol1(xdep, errarr2, effdep2, ijob, inter_type) ;
602 for (
int i=0; i<nxline; i++) {
603 while (fabs(xlk[i].
x - xdep(iresu)) > x12 ) {
608 farr.
set(it,ir) = effdep2(iresu) ;
613 for (
int i=0; i<nxline; i++) {
614 resu = effarr2(iresu) ;
616 if ((xlk[i+1].
x-xlk[i].
x) > x12) {
622 farr.
set(it,ir) = resu ;
639 const Tbl& parr,
const int type_inter)
const {
653 Tbl *fdept =
new Tbl(ntv, nrv) ;
655 Tbl intermediaire(npv, ntm, nrm) ;
658 Tbl farr(npm, ntm, nrm) ;
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] ;
665 for (
int j=0; j<ntm; j++)
666 for (
int k=0; k<nrm; k++) intermediaire.
set(i,j,k) = fr(j,k) ;
671 fdept =
new Tbl(npv) ;
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) ;
678 for (
int k=0; k<npm; k++) farr.
set(k,i,j) = fr(k) ;
686 const Tbl& tarr,
const Tbl& parr,
const
687 int inter_type)
const {
700 Tbl farr(np, nt, nr) ;
703 bool coq = (rarr(0)/rarr(nr-1) > 1.e-6) ;
706 rarr2 =
new Tbl(2*nr) ;
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) ;
713 int nr2 = coq ? 2*nr : nr ;
715 Tbl cylindre(nz, np, nr2) ;
717 for(
int iz=0; iz<nz; iz++) {
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) ;
724 cercle =
interpol2c(*
x, *
y, carre, coq ? *rarr2 : rarr, parr, inter_type) ;
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) ;
731 for (
int ip=0; ip<np; ip++) {
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,
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) ;
745 if (coq)
delete rarr2 ;
Active physical coordinates and mapping derivatives.
int * dim
Array of dimensions (size: ndim).
int ndim
Number of dimensions of the Tbl: can be 1, 2 or 3.
int type_p
Type of symmetry in :
int nfantome
The number of hidden cells (same on each side)
double * zrmin
Lower boundary for z (or r ) direction
Dim_tbl dim
The dimensions of the grid.
Tbl * zr
Arrays containing the values of coordinate z (or r) on the nodes
double get_zr(const int i) const
Read-only of a particular value of the coordinate z (or r ) at the nodes.
int type_t
Type of symmetry in :
Tbl interpol1(const Tbl &rdep, const Tbl &rarr, const Tbl &fdep, int flag, const int type_inter) const
Performs 1D interpolation.
double * zrmax
Higher boundary for z (or r ) direction
double * xmax
Higher boundary for x dimension.
double * xmin
Lower boundary for x dimension.
double * ymin
Lower boundary for y dimension.
double * ymax
Higher boundary for y dimension.
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).
virtual Tbl interpol3(const Tbl &fdep, const Tbl &rarr, const Tbl &tetarr, const Tbl &phiarr, const int type_inter) const
Performs 3D interpolation.
Tbl * y
Arrays containing the values of coordinate y on the nodes.
virtual Tbl interpol2(const Tbl &fdep, const Tbl &rarr, const Tbl &tetarr, const int type_inter) const
Performs 2D interpolation.
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...
Tbl * x
Arrays containing the values of coordinate x on the nodes.
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,...
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...
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,...
Tbl * phi
Arrays containing the values of coordinate on the nodes.
virtual Tbl interpol3(const Tbl &fdep, const Tbl &rarr, const Tbl &tetarr, const Tbl &phiarr, const int type_inter) const
Performs 3D interpolation.
virtual Tbl interpol2(const Tbl &fdep, const Tbl &rarr, const Tbl &tetarr, const int type_inter) const
Performs 2D interpolation.
Tbl * tet
Arrays containing the values of coordinate on the nodes.
Base class for coordinate mappings.
Coord r
r coordinate centered on the grid
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Dim_tbl dim
Number of dimensions, size,...
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double & set(int i)
Read/write of a particular element (index i) (1D case)
double * t
The array of double.
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Cmp sin(const Cmp &)
Sine.
Cmp cos(const Cmp &)
Cosine.