LORENE
map_af_primr.C
1 /*
2  * Method Map_af::primr
3  *
4  * (see file map.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2004 Eric Gourgoulhon
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 version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 char map_af_primr_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_af_primr.C,v 1.7 2014/10/13 08:53:03 j_novak Exp $" ;
29 
30 /*
31  * $Id: map_af_primr.C,v 1.7 2014/10/13 08:53:03 j_novak Exp $
32  * $Log: map_af_primr.C,v $
33  * Revision 1.7 2014/10/13 08:53:03 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.6 2014/10/06 15:13:12 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.5 2013/04/25 15:46:05 j_novak
40  * Added special treatment in the case np = 1, for type_p = NONSYM.
41  *
42  * Revision 1.4 2007/12/20 09:11:05 jl_cornou
43  * Correction of an error in op_sxpun about Jacobi(0,2) polynomials
44  *
45  * Revision 1.3 2004/07/26 16:02:23 j_novak
46  * Added a flag to specify whether the primitive should be zero either at r=0
47  * or at r going to infinity.
48  *
49  * Revision 1.1 2004/06/14 15:25:34 e_gourgoulhon
50  * First version.
51  *
52  *
53  * $Header: /cvsroot/Lorene/C++/Source/Map/map_af_primr.C,v 1.7 2014/10/13 08:53:03 j_novak Exp $
54  *
55  */
56 
57 
58 // C headers
59 #include <cstdlib>
60 
61 // Lorene headers
62 #include "map.h"
63 #include "tensor.h"
64 
65 namespace Lorene {
66 void _primr_pas_prevu(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
67 void _primr_r_cheb(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
68 void _primr_r_chebp(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
69 void _primr_r_chebi(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
70 void _primr_r_chebpim_p(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
71 void _primr_r_chebpim_i(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
72 void _primr_r_jaco02(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
73 
74 void Map_af::primr(const Scalar& uu, Scalar& resu, bool null_infty) const {
75 
76  static void (*prim_domain[MAX_BASE])(const Tbl&, int bin, const Tbl&,
77  Tbl&, int&, Tbl& ) ;
78  static bool first_call = true ;
79 
80  // Initialisation at first call of the array of primitive functions
81  // depending upon the basis in r
82  // ----------------------------------------------------------------
83  if (first_call) {
84  for (int i=0 ; i<MAX_BASE ; i++) prim_domain[i] = _primr_pas_prevu ;
85 
86  prim_domain[R_CHEB >> TRA_R] = _primr_r_cheb ;
87  prim_domain[R_CHEBU >> TRA_R] = _primr_r_cheb ;
88  prim_domain[R_CHEBP >> TRA_R] = _primr_r_chebp ;
89  prim_domain[R_CHEBI >> TRA_R] = _primr_r_chebi ;
90  prim_domain[R_CHEBPIM_P >> TRA_R] = _primr_r_chebpim_p ;
91  prim_domain[R_CHEBPIM_I >> TRA_R] = _primr_r_chebpim_i ;
92  prim_domain[R_JACO02 >> TRA_R] = _primr_r_jaco02 ;
93 
94  first_call = false ;
95  }
96 
97  // End of first call operations
98  // ----------------------------
99 
100  assert(uu.get_etat() != ETATNONDEF) ;
101  assert(uu.get_mp().get_mg() == mg) ;
102  assert(resu.get_mp().get_mg() == mg) ;
103 
104  // Special case of vanishing input:
105  if (uu.get_etat() == ETATZERO) {
106  resu.set_etat_zero() ;
107  return ;
108  }
109 
110  // General case
111  assert( (uu.get_etat() == ETATQCQ) || (uu.get_etat() == ETATUN) ) ;
112  assert(uu.check_dzpuis(2)) ;
113 
114  int nz = mg->get_nzone() ;
115  int nzm1 = nz - 1 ;
116  int np = mg->get_np(0) ;
117  int nt = mg->get_nt(0) ;
118 #ifndef NDEBUG
119  for (int l=1; l<nz; l++) {
120  assert (mg->get_np(l) == np) ;
121  assert (mg->get_nt(l) == nt) ;
122  }
123 #endif
124 
125  const Valeur& vuu = uu.get_spectral_va() ;
126  vuu.coef() ;
127 
128  const Mtbl_cf& cuu = *(vuu.c_cf) ;
129  assert(cuu.t != 0x0) ;
130 
131  const Base_val& buu = vuu.get_base() ; // spectral bases of the input
132 
133  resu.set_etat_qcq() ; // result in ordinary state
134  Valeur& vprim = resu.set_spectral_va() ;
135  vprim.set_etat_cf_qcq() ; // allocates the Mtbl_cf for the coefficients
136  // of the result
137  Mtbl_cf& cprim = *(vprim.c_cf) ;
138  cprim.set_etat_qcq() ; // allocates the Tbl's to store the coefficients
139  // of the result in each domain
140 
141  Base_val& bprim = cprim.base ; // spectral bases of the result
142 
143  Tbl val_rmin(np+2,nt) ; // Values of primitive at the left boundary
144  // in the current domain
145  Tbl val_rmax(np+2,nt) ; // same but for the right boundary
146 
147  val_rmin.set_etat_zero() ; // initialisation: primitive = 0 at r=0
148 
149  int lmax = (mg->get_type_r(nzm1) == UNSURR) ? nz-2 : nzm1 ;
150 
151  for (int l=0; l<=lmax; l++) {
152  assert(cuu.t[l] != 0x0) ;
153  assert(cprim.t[l] != 0x0) ;
154  const Tbl& cfuu = *(cuu.t[l]) ;
155  Tbl& cfprim = *(cprim.t[l]) ;
156 
157  int buu_dom = buu.get_b(l) ;
158  int base_r = (buu_dom & MSQ_R) >> TRA_R ;
159 
160  prim_domain[base_r](cfuu, buu_dom, val_rmin, cfprim, bprim.b[l],
161  val_rmax) ;
162 
163  cfprim *= alpha[l] ;
164  val_rmin = alpha[l] * val_rmax / alpha[l+1] ; // for next domain
165  }
166 
167  // Special case of compactified external domain (CED)
168  // --------------------------------------------------
169  if (mg->get_type_r(nzm1) == UNSURR) {
170  val_rmin = - val_rmin ;
171  const Tbl& cfuu = *(cuu.t[nzm1]) ;
172  Tbl& cfprim = *(cprim.t[nzm1]) ;
173 
174  int buu_dom = buu.get_b(nzm1) ;
175  int base_r = (buu_dom & MSQ_R) >> TRA_R ;
176  assert(base_r == R_CHEBU) ;
177 
178  prim_domain[base_r](cfuu, buu_dom, val_rmin, cfprim, bprim.b[nzm1],
179  val_rmax) ;
180 
181  cfprim *= - alpha[nzm1] ;
182  }
183 
184  if (null_infty)
185  for (int k=0; k<np; k++) //## not very elegant!
186  for(int j=0; j<nt; j++)
187  val_rmax.set(k,j) = cprim.val_out_bound_jk(nzm1, j, k) ;
188 
189  // The output spectral bases (set on the Mtbl_cf) are copied to the Valeur:
190  vprim.set_base(bprim) ;
191 
192  if (null_infty)
193  for (int l=0; l<nz; l++) //## not very elegant!
194  for (int k=0; k<np; k++)
195  for(int j=0; j<nt; j++)
196  for (int i=0; i<mg->get_nr(l); i++)
197  vprim.set(l, k, j, i) -= val_rmax(k,j) ;
198 
199 
200 }
201 }
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
virtual void primr(const Scalar &uu, Scalar &resu, bool null_infty) const
Computes the radial primitive which vanishes for .
Definition: map_af_primr.C:74
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
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 coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:324
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition: map.h:2033
Lorene prototypes.
Definition: app_hor.h:64
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:353
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition: valeur.C:810
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
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...
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
int * b
Array (size: nzone ) of the spectral basis in each domain.
Definition: base_val.h:331
#define MSQ_R
Extraction de l&#39;info sur R.
Definition: type_parite.h:152
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl_cf.C:300
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
Definition: type_parite.h:178
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
Definition: type_parite.h:176
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: tbl.C:347
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined.
Definition: map.h:676
Bases of the spectral expansions.
Definition: base_val.h:322
int get_b(int l) const
Returns the code for the expansion basis in domain no. l.
Definition: base_val.h:389
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:186
Base_val base
Bases of the spectral expansions.
Definition: mtbl_cf.h:200
Basic array class.
Definition: tbl.h:161
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:474
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition: scalar.C:873
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:861
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s which contain the spectral coefficients in each domain...
Definition: mtbl_cf.h:205
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:601
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition: valeur.h:363
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166