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
30char 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
93namespace Lorene {
94Scalar 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 ;
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 }
196 Tbl rarr(nrarr) ;
197 rarr.set_etat_qcq() ;
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 ;
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 }
249 Tbl rarr(nrarr) ;
250 rarr.set_etat_qcq() ;
252 tetarr.set_etat_qcq() ;
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
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}
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Class for spherical Godunov-type grids.
Definition grille_val.h:579
Base class for coordinate mappings.
Definition map.h:670
Coord r
r coordinate centered on the grid
Definition map.h:718
Coord tet
coordinate centered on the grid
Definition map.h:719
Coord phi
coordinate centered on the grid
Definition map.h:720
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
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Grille_val * gval
The Grille_val (cartesian or spherical) on which the array is defined.
Definition tbl_val.h:110
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).
const Dim_tbl * dim
The Dim_tbl giving the dimensions and number of points (without the hidden cells).
Definition tbl_val.h:108
double * txti
The array at x (or ) interfaces.
Definition tbl_val.h:118
double * tzri
The array at z (or r) interfaces.
Definition tbl_val.h:116
void annule_hard()
Sets the Tbl_val to zero in a hard way.
Definition tbl_val.C:311
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl_val.C:294
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 .
double * t
The array of double at the nodes.
Definition tbl_val.h:114
int get_taille() const
Gives the size of the node array (including the hidden cells)
Definition tbl_val.h:462
int etat
logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition tbl_val.h:103
Basic array class.
Definition tbl.h:161
Lorene prototypes.
Definition app_hor.h:64