LORENE
scalar_manip.C
1/*
2 * Copyright (c) 2003 Eric Gourgoulhon & Jerome Novak
3 *
4 * Copyright (c) 2000-2001 Philippe Grandclement (Cmp version)
5 *
6 * This file is part of LORENE.
7 *
8 * LORENE is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * LORENE is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with LORENE; if not, write to the Free Software
20 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 *
22 */
23
24
25char scalar_manip_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_manip.C,v 1.19 2014/10/13 08:53:46 j_novak Exp $" ;
26
27/*
28 * $Id: scalar_manip.C,v 1.19 2014/10/13 08:53:46 j_novak Exp $
29 * $Log: scalar_manip.C,v $
30 * Revision 1.19 2014/10/13 08:53:46 j_novak
31 * Lorene classes and functions now belong to the namespace Lorene.
32 *
33 * Revision 1.18 2014/10/06 15:16:15 j_novak
34 * Modified #include directives to use c++ syntax.
35 *
36 * Revision 1.17 2008/10/03 09:03:52 j_novak
37 * Correction of yet another mistake (the array values in physical space was not
38 * destroyed).
39 *
40 * Revision 1.16 2008/09/30 08:35:18 j_novak
41 * Correction of forgotten call to coef()
42 *
43 * Revision 1.15 2008/09/29 13:23:51 j_novak
44 * Implementation of the angular mapping associated with an affine
45 * mapping. Things must be improved to take into account the domain index.
46 *
47 * Revision 1.14 2008/09/22 19:08:01 j_novak
48 * New methods to deal with boundary conditions
49 *
50 * Revision 1.13 2007/06/21 20:00:00 k_taniguchi
51 * Addition of the method filtre_r (int n, int nz)
52 *
53 * Revision 1.12 2006/06/28 07:46:41 j_novak
54 * Better treatment in the case of a domain set to zero.
55 *
56 * Revision 1.11 2005/09/07 13:39:10 j_novak
57 * *** empty log message ***
58 *
59 * Revision 1.10 2005/09/07 13:10:48 j_novak
60 * Added a filter setting to zero all mulitpoles in a given range.
61 *
62 * Revision 1.9 2004/11/23 12:47:44 f_limousin
63 * Add function filtre_tp(int nn, int nz1, int nz2).
64 *
65 * Revision 1.8 2004/05/07 11:24:54 f_limousin
66 * Implement new method filtre_r (int* nn)
67 *
68 * Revision 1.7 2004/02/27 09:47:26 f_limousin
69 * New methods filtre_phi(int) and filtre_theta(int).
70 *
71 * Revision 1.6 2004/01/23 13:26:28 e_gourgoulhon
72 * Added methods set_inner_boundary and set_outer_boundary.
73 * Methods set_val_inf and set_val_hor, which are particular cases of
74 * the above, have been suppressed.
75 *
76 * Revision 1.5 2003/11/13 13:43:55 p_grandclement
77 * Addition of things needed for Bhole::update_metric (const Etoile_bin&, double, double)
78 *
79 * Revision 1.4 2003/10/11 14:47:17 e_gourgoulhon
80 * Lines 112 and 145 : replaced "0" by "double(0)" in the comparison
81 * statement.
82 *
83 * Revision 1.3 2003/10/10 15:57:29 j_novak
84 * Added the state one (ETATUN) to the class Scalar
85 *
86 * Revision 1.2 2003/10/08 14:24:09 j_novak
87 * replaced mult_r_zec with mult_r_ced
88 *
89 * Revision 1.1 2003/09/25 09:33:36 j_novak
90 * Added methods for integral calculation and various manipulations
91 *
92 *
93 * $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_manip.C,v 1.19 2014/10/13 08:53:46 j_novak Exp $
94 *
95 */
96
97//standard
98#include <cstdlib>
99#include <cmath>
100
101// Lorene
102#include "tensor.h"
103#include "proto.h"
104#include "utilitaires.h"
105
106/*
107 * Annule tous les l entre l_min et l_max (compris)
108 */
109
110namespace Lorene {
111void Scalar::annule_l (int l_min, int l_max, bool ylm_output) {
112
113 assert (etat != ETATNONDEF) ;
114 assert (l_min <= l_max) ;
115 assert (l_min >= 0) ;
116 if (etat == ETATZERO )
117 return ;
118
119 if (etat == ETATUN) {
120 if (l_min == 0) set_etat_zero() ;
121 else return ;
122 }
123
124 va.ylm() ;
125 Mtbl_cf& m_coef = *va.c_cf ;
126 const Base_val& base = va.base ;
127 int l_q, m_q, base_r ;
128 for (int lz=0; lz<mp->get_mg()->get_nzone(); lz++)
129 if (m_coef(lz).get_etat() != ETATZERO)
130 for (int k=0; k<mp->get_mg()->get_np(lz)+1; k++)
131 for (int j=0; j<mp->get_mg()->get_nt(lz); j++)
132 for (int i=0; i<mp->get_mg()->get_nr(lz); i++) {
133 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
134 if ((l_min <= l_q) && (l_q<= l_max))
135 m_coef.set(lz, k, j, i) = 0 ;
136 }
137 if (va.c != 0x0) {
138 delete va.c ;
139 va.c = 0x0 ;
140 }
141 if (!ylm_output)
142 va.ylm_i() ;
143
144 return ;
145}
146
147/*
148 * Annule les n derniers coefficients en r dans la derniere zone
149 */
150
151void Scalar::filtre (int n) {
152
153 assert (etat != ETATNONDEF) ;
154 if ( (etat == ETATZERO) || (etat == ETATUN) )
155 return ;
156
157 int nz = mp->get_mg()->get_nzone() ;
158 int np = mp->get_mg()->get_np(nz-1) ;
159 int nt = mp->get_mg()->get_nt(nz-1) ;
160 int nr = mp->get_mg()->get_nr(nz-1) ;
161
162 del_deriv() ;
163
164 va.coef() ;
166
167
168 for (int k=0 ; k<np+1 ; k++)
169 if (k!=1)
170 for (int j=0 ; j<nt ; j++)
171 for (int i=nr-1 ; i>nr-1-n ; i--)
172 va.c_cf->set(nz-1, k, j, i) = 0 ;
173}
174
175
176/*
177 * Annule les n derniers coefficients en r dans toutes les zones
178 */
179
180void Scalar::filtre_r (int* nn) {
181 assert (etat != ETATNONDEF) ;
182 if ( (etat == ETATZERO) || (etat == ETATUN) )
183 return ;
184
185 del_deriv() ;
186
187 va.coef() ;
189 int nz = mp->get_mg()->get_nzone() ;
190 int* nr = new int[nz];
191 int* nt = new int[nz];
192 int* np = new int[nz];
193 for (int l=0; l<=nz-1; l++) {
194 nr[l] = mp->get_mg()->get_nr(l) ;
195 nt[l] = mp->get_mg()->get_nt(l) ;
196 np[l] = mp->get_mg()->get_np(l) ;
197 }
198
199 for (int l=0; l<=nz-1; l++)
200 for (int k=0 ; k<np[l]+1 ; k++)
201 if (k!=1)
202 for (int j=0 ; j<nt[l] ; j++)
203 for (int i=nr[l]-1; i>nr[l]-1-nn[l] ; i--)
204 va.c_cf->set(l, k, j, i) = 0 ;
205
206 if (va.c != 0x0) {
207 delete va.c ;
208 va.c = 0x0 ;
209 }
210
211}
212
213
214/*
215 * Annule les n derniers coefficients en r dans zone nz
216 */
217
218void Scalar::filtre_r (int n, int nz) {
219 assert (etat != ETATNONDEF) ;
220 if ( (etat == ETATZERO) || (etat == ETATUN) )
221 return ;
222
223 del_deriv() ;
224
225 va.coef() ;
227 int nr = mp->get_mg()->get_nr(nz) ;
228 int nt = mp->get_mg()->get_nt(nz) ;
229 int np = mp->get_mg()->get_np(nz) ;
230
231 for (int k=0 ; k<np+1 ; k++)
232 if (k!=1)
233 for (int j=0 ; j<nt ; j++)
234 for (int i=nr-1; i>nr-1-n ; i--)
235 va.c_cf->set(nz, k, j, i) = 0 ;
236
237 if (va.c != 0x0) {
238 delete va.c ;
239 va.c = 0x0 ;
240 }
241
242}
243
244
245/*
246 * Annule les n derniers coefficients en phi dans zone nz
247 */
248
249void Scalar::filtre_phi (int n, int nz) {
250 assert (etat != ETATNONDEF) ;
251 if ( (etat == ETATZERO) || (etat == ETATUN) )
252 return ;
253
254 del_deriv() ;
255
256 va.coef() ;
258 int np = mp->get_mg()->get_np(nz) ;
259 int nt = mp->get_mg()->get_nt(nz) ;
260 int nr = mp->get_mg()->get_nr(nz) ;
261
262 for (int k=np+1-n ; k<np+1 ; k++)
263 for (int j=0 ; j<nt ; j++)
264 for (int i=0 ; i<nr ; i++)
265 va.c_cf->set(nz, k, j, i) = 0 ;
266
267}
268
269
270void Scalar::filtre_tp(int nn, int nz1, int nz2) {
271
272 va.filtre_tp(nn, nz1, nz2) ;
273
274}
275
276
277 /* Sets the value of the {\tt Scalar} at the inner boundary of a given
278 * domain.
279 * @param l [input] domain index
280 * @param x [input] (constant) value at the inner boundary of domain no. {\tt l}
281 */
282
284
285 assert (etat != ETATNONDEF) ;
286 if (etat == ETATZERO) {
287 if (x0 == double(0)) return ;
288 else annule_hard() ;
289 }
290
291 if (etat == ETATUN) {
292 if (x0 == double(1)) return ;
293 else etat = ETATQCQ ;
294 }
295
296 del_deriv() ;
297
298 int nt = mp->get_mg()->get_nt(l0) ;
299 int np = mp->get_mg()->get_np(l0) ;
300
301 va.coef_i() ;
303
304 for (int k=0 ; k<np ; k++)
305 for (int j=0 ; j<nt ; j++)
306 va.set(l0, k, j, 0) = x0 ;
307}
308
309 /* Sets the value of the {\tt Scalar} at the outer boundary of a given
310 * domain.
311 * @param l [input] domain index
312 * @param x [input] (constant) value at the outer boundary of domain no. {\tt l}
313 */
314
316
317 assert (etat != ETATNONDEF) ;
318 if (etat == ETATZERO) {
319 if (x0 == double(0)) return ;
320 else annule_hard() ;
321 }
322
323 if (etat == ETATUN) {
324 if (x0 == double(1)) return ;
325 else etat = ETATQCQ ;
326 }
327
328 del_deriv() ;
329
330 int nrm1 = mp->get_mg()->get_nr(l0) - 1 ;
331 int nt = mp->get_mg()->get_nt(l0) ;
332 int np = mp->get_mg()->get_np(l0) ;
333
334 va.coef_i() ;
336
337 for (int k=0 ; k<np ; k++)
338 for (int j=0 ; j<nt ; j++)
339 va.set(l0, k, j, nrm1) = x0 ;
340}
341
342/*
343 * Permet de fixer la decroissance du cmp a l infini en viurant les
344 * termes en 1/r^n
345 */
347
348 if (puis<dzpuis)
349 return ;
350 else {
351
352 int nbre = puis-dzpuis ;
353
354 // le confort avant tout ! (c'est bien le confort ...)
355 int nz = mp->get_mg()->get_nzone() ;
356 int np = mp->get_mg()->get_np(nz-1) ;
357 int nt = mp->get_mg()->get_nt(nz-1) ;
358 int nr = mp->get_mg()->get_nr(nz-1) ;
359
360 const Map_af* map = dynamic_cast<const Map_af*>(mp) ;
361 if (map == 0x0) {
362 cout << "Le mapping doit etre affine" << endl ;
363 abort() ;
364 }
365
366 double alpha = map->get_alpha()[nz-1] ;
367
368 Scalar courant (*this) ;
369
370 va.coef() ;
372
373 for (int conte=0 ; conte<nbre ; conte++) {
374
375 int base_r = courant.va.base.get_base_r(nz-1) ;
376
377 courant.va.coef() ;
378
379 // On calcul les coefficients de 1/r^conte
380 double* coloc = new double [nr] ;
381 int * deg = new int[3] ;
382 deg[0] = 1 ;
383 deg[1] = 1 ;
384 deg[2] = nr ;
385
386 for (int i=0 ; i<nr ; i++)
387 coloc[i] =pow(alpha, double(conte))*
388 pow(-1-cos(M_PI*i/(nr-1)), double(conte)) ;
389
390 cfrcheb(deg, deg, coloc, deg, coloc) ;
391
392 for (int k=0 ; k<np+1 ; k++)
393 if (k != 1)
394 for (int j=0 ; j<nt ; j++) {
395
396 // On doit determiner le coefficient du truc courant :
397 double* coef = new double [nr] ;
398 double* auxi = new double[1] ;
399 for (int i=0 ; i<nr ; i++)
400 coef[i] = (*courant.va.c_cf)(nz-1, k, j, i) ;
401 switch (base_r) {
402 case R_CHEBU :
403 som_r_chebu (coef, nr, 1, 1, 1, auxi) ;
404 break ;
405 default :
406 som_r_pas_prevu (coef, nr, 1, 1, 1, auxi) ;
407 break ;
408 }
409
410 // On modifie le cmp courant :
411 courant.va.coef() ;
412 courant.va.set_etat_cf_qcq() ;
413 courant.va.c_cf->set(nz-1, k, j, 0) -= *auxi ;
414
415 for (int i=0 ; i<nr ; i++)
416 this->va.c_cf->set(nz-1, k, j, i) -= *auxi * coloc[i] ;
417
418
419 delete [] coef ;
420 delete [] auxi ;
421 }
422 delete [] coloc ;
423 delete [] deg ;
424
425 courant.mult_r_ced() ;
426 }
427 }
428}
429
431 va.coef() ;
432 if (output_ylm) va.ylm() ;
433
434 int np = mp->get_mg()->get_np(l_zone) ;
435 int nt = mp->get_mg()->get_nt(l_zone) ;
436 Tbl resu(np+2, nt) ;
437 if (etat == ETATZERO) resu.set_etat_zero() ;
438 else {
439 assert(etat == ETATQCQ) ;
440 resu.set_etat_qcq() ;
441 for (int k=0; k<np+2; k++)
442 for (int j=0; j<nt; j++)
443 resu.set(k, j) = va.c_cf->val_out_bound_jk(l_zone, j, k) ;
444 }
445 return resu ;
446}
447
449 assert(mp->get_mg()->get_type_r(l_zone) != RARE) ;
450 va.coef() ;
451 if (output_ylm) va.ylm() ;
452
453 int np = mp->get_mg()->get_np(l_zone) ;
454 int nt = mp->get_mg()->get_nt(l_zone) ;
455 Tbl resu(np+2, nt) ;
456 if (etat == ETATZERO) resu.set_etat_zero() ;
457 else {
458 assert(etat == ETATQCQ) ;
459 resu.set_etat_qcq() ;
460 for (int k=0; k<np+2; k++)
461 for (int j=0; j<nt; j++)
462 resu.set(k, j) = va.c_cf->val_in_bound_jk(l_zone, j, k) ;
463 }
464 return resu ;
465}
466
468 va.coef() ;
469 if (output_ylm) va.ylm() ;
470
471 Scalar resu(mp->mp_angu(l_zone)) ;
472 resu.std_spectral_base() ;
473 Base_val base = resu.get_spectral_base() ;
475 resu.set_spectral_base(base) ;
476 if (etat == ETATZERO) resu.set_etat_zero() ;
477 else {
478 assert(etat == ETATQCQ) ;
479 resu.annule_hard() ;
480 int np = mp->get_mg()->get_np(l_zone) ;
481 int nt = mp->get_mg()->get_nt(l_zone) ;
482 for (int k=0; k<np+2; k++)
483 for (int j=0; j<nt; j++)
484 resu.set_spectral_va().c_cf->set(0, k, j, 0)
486 delete resu.set_spectral_va().c ;
487 resu.set_spectral_va().c = 0x0 ;
488 }
489 return resu ;
490}
491}
Bases of the spectral expansions.
Definition base_val.h:322
void set_base_t(int base_t)
Sets the expansion basis for functions in all domains.
Definition base_val.C:195
int get_base_t(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
Definition base_val.h:411
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Affine radial mapping.
Definition map.h:2027
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition map_af.C:477
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:294
double val_out_bound_jk(int l, int j, int k) const
Computes the angular coefficient of index j,k of the field represented by *this at by means of the s...
double val_in_bound_jk(int l, int j, int k) const
Computes the angular coefficient of index j,k of the field represented by *this at by means of the s...
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
void annule_l(int l_min, int l_max, bool ylm_output=false)
Sets all the multipolar components between l_min and l_max to zero.
virtual void del_deriv() const
Logical destructor of the derivatives.
Definition scalar.C:287
friend Scalar cos(const Scalar &)
Cosine.
Scalar scalar_out_bound(int n, bool leave_ylm=false)
Returns the Scalar containing the values of angular coefficients at the outer boundary.
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
void fixe_decroissance(int puis)
Substracts all the components behaving like in the external domain, with n strictly lower than puis ...
Tbl tbl_in_bound(int n, bool leave_ylm=false)
Returns the Tbl containing the values of angular coefficients at the inner boundary.
Tbl tbl_out_bound(int l_dom, bool leave_ylm=false)
Returns the Tbl containing the values of angular coefficients at the outer boundary.
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
void filtre_tp(int nn, int nz1, int nz2)
Sets the n lasts coefficients in to 0 in the domain nz1 to nz2 when expressed in spherical harmonics...
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary).
Definition scalar.h:396
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
void set_inner_boundary(int l, double x)
Sets the value of the Scalar at the inner boundary of a given domain.
void filtre_r(int *nn)
Sets the n lasts coefficients in r to 0 in all domains.
friend Scalar pow(const Scalar &, int)
Power .
void filtre_phi(int n, int zone)
Sets the n lasts coefficients in to 0 in the domain zone .
void set_outer_boundary(int l, double x)
Sets the value of the Scalar at the outer boundary of a given domain.
Valeur va
The numerical value of the Scalar
Definition scalar.h:405
int dzpuis
Power of r by which the quantity represented by this must be divided in the compactified external d...
Definition scalar.h:403
Basic array class.
Definition tbl.h:161
void filtre_tp(int nn, int nz1, int nz2)
Sets the n lasts coefficients in to 0 in the domain nz1 to nz2 when expressed in spherical harmonics...
Definition valeur.C:924
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
Definition valeur.C:712
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c ).
Definition valeur.C:701
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Mtbl * c
Values of the function at the points of the multi-grid
Definition valeur.h:299
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition valeur.h:363
void coef_i() const
Computes the physical value of *this.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
void coef() const
Computes the coeffcients of *this.
void ylm_i()
Inverse of ylm()
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Lorene prototypes.
Definition app_hor.h:64