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
28char 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
65namespace Lorene {
66void _primr_pas_prevu(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
67void _primr_r_cheb(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
68void _primr_r_chebp(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
69void _primr_r_chebi(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
70void _primr_r_chebpim_p(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
71void _primr_r_chebpim_i(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
72void _primr_r_jaco02(const Tbl&, int, const Tbl&, Tbl&, int&, Tbl&) ;
73
74void 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}
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
int * b
Array (size: nzone ) of the spectral basis in each domain.
Definition base_val.h:331
virtual void primr(const Scalar &uu, Scalar &resu, bool null_infty) const
Computes the radial primitive which vanishes for .
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition map.h:2033
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition map.h:676
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:474
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
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
Definition mtbl_cf.h:205
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition mtbl_cf.C:300
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...
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
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
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition scalar.h:601
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
Basic array class.
Definition tbl.h:161
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition tbl.C:347
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
Values and coefficients of a (real-value) function.
Definition valeur.h:287
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
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Definition valeur.h:480
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition valeur.C:810
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition valeur.h:363
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
void coef() const
Computes the coeffcients of *this.
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define MSQ_R
Extraction de l'info sur R.
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
#define R_CHEB
base de Chebychev ordinaire (fin)
#define R_CHEBP
base de Cheb. paire (rare) seulement
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene prototypes.
Definition app_hor.h:64