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
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 $" ;
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
96namespace Lorene {
97bool 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
168bool 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 )
214bool 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 )
274bool 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
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) ;
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
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() ;
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] ;
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] ;
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
456struct Point {
457 double x ;
458 int l ;
459 int k ;
460};
461
462#endif /* DOXYGEN_SHOULD_SKIP_THIS */
463
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 ;
468}
469
471 const Tbl& tetarr, const int type_inter) const
472{
473 return interpol2c(*zr, *x, fdep, rarr, tetarr, type_inter) ;
474}
475
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) ;
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
581 errarr2.set_etat_qcq() ;
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 ;
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
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() ;
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] ;
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) ;
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
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...
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) ;
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}
Active physical coordinates and mapping derivatives.
Definition coord.h:90
int * dim
Array of dimensions (size: ndim).
Definition dim_tbl.h:102
int ndim
Number of dimensions of the Tbl: can be 1, 2 or 3.
Definition dim_tbl.h:101
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
int type_p
Type of symmetry in :
Definition grille_val.h:114
int nfantome
The number of hidden cells (same on each side)
Definition grille_val.h:104
double * zrmin
Lower boundary for z (or r ) direction
Definition grille_val.h:117
Dim_tbl dim
The dimensions of the grid.
Definition grille_val.h:102
Tbl * zr
Arrays containing the values of coordinate z (or r) on the nodes
Definition grille_val.h:124
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
int type_t
Type of symmetry in :
Definition grille_val.h:109
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
Definition grille_val.h:120
double * xmax
Higher boundary for x dimension.
Definition grille_val.h:335
double * xmin
Lower boundary for x dimension.
Definition grille_val.h:333
double * ymin
Lower boundary for y dimension.
Definition grille_val.h:337
double * ymax
Higher boundary for y dimension.
Definition grille_val.h:339
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.
Definition grille_val.h:347
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.
Definition grille_val.h:343
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.
Definition grille_val.h:591
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.
Definition grille_val.h:587
Affine radial mapping.
Definition map.h:2027
Base class for coordinate mappings.
Definition map.h:670
Coord r
r coordinate centered on the grid
Definition map.h:718
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Multi-domain grid.
Definition grilles.h:273
Basic array class.
Definition tbl.h:161
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition tbl.h:403
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:69
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
Lorene prototypes.
Definition app_hor.h:64