LORENE
map_af_integ.C
1/*
2 * Method of the class Map_af for computing the integral of a Cmp over
3 * all space.
4 */
5
6/*
7 * Copyright (c) 1999-2003 Eric Gourgoulhon
8 *
9 * This file is part of LORENE.
10 *
11 * LORENE is free software; you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License as published by
13 * the Free Software Foundation; either version 2 of the License, or
14 * (at your option) any later version.
15 *
16 * LORENE is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with LORENE; if not, write to the Free Software
23 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24 *
25 */
26
27
28char map_af_integ_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_af_integ.C,v 1.9 2014/10/13 08:53:02 j_novak Exp $" ;
29
30/*
31 * $Id: map_af_integ.C,v 1.9 2014/10/13 08:53:02 j_novak Exp $
32 * $Log: map_af_integ.C,v $
33 * Revision 1.9 2014/10/13 08:53:02 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.8 2014/10/06 15:13:12 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.7 2009/10/08 16:20:47 j_novak
40 * Addition of new bases T_COS and T_SIN.
41 *
42 * Revision 1.6 2008/08/27 08:48:05 jl_cornou
43 * Added R_JACO02 case
44 *
45 * Revision 1.5 2007/10/05 15:56:19 j_novak
46 * Addition of new bases for the non-symmetric case in theta.
47 *
48 * Revision 1.4 2003/12/19 16:21:43 j_novak
49 * Shadow hunt
50 *
51 * Revision 1.3 2003/10/19 19:58:15 e_gourgoulhon
52 * Access to Base_val::b now via Base_val::get_b().
53 *
54 * Revision 1.2 2003/10/15 10:35:27 e_gourgoulhon
55 * Changed cast (double *) into static_cast<double*>.
56 *
57 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
58 * LORENE
59 *
60 * Revision 1.2 2000/01/28 16:09:37 eric
61 * Remplacement du ci.get_dzpuis() == 4 par ci.check_dzpuis(4).
62 *
63 * Revision 1.1 1999/12/09 10:45:43 eric
64 * Initial revision
65 *
66 *
67 * $Header: /cvsroot/Lorene/C++/Source/Map/map_af_integ.C,v 1.9 2014/10/13 08:53:02 j_novak Exp $
68 *
69 */
70
71// Headers C
72#include <cstdlib>
73#include <cmath>
74
75
76// Headers Lorene
77#include "map.h"
78#include "cmp.h"
79
80namespace Lorene {
81Tbl* Map_af::integrale(const Cmp& ci) const {
82
83 static double* cx_tcp = 0 ; // Coefficients theta, dev. en cos(2l theta)
84 static double* cx_rrp = 0 ; // Coefficients r rare, dev. en T_{2i}
85 static double* cx_rf_x2 = 0 ; // Coefficients r fin, int. en r^2
86 static double* cx_rf_x = 0 ; // Coefficients r fin, int en r
87 static double* cx_rf = 0 ; // Coefficients r fin, int en cst.
88
89 static int nt_cp_pre = 0 ;
90 static int nr_p_pre = 0 ;
91 static int nr_f_pre = 0 ;
92
93 assert(ci.get_etat() != ETATNONDEF) ;
94
95 int nz = mg->get_nzone() ;
96
97 Tbl* resu = new Tbl(nz) ;
98
99 if (ci.get_etat() == ETATZERO) {
100 resu->annule_hard() ;
101 return resu ;
102 }
103
104 assert( ci.get_etat() == ETATQCQ ) ;
105
106 (ci.va).coef() ; // The spectral coefficients are required
107
108 const Mtbl_cf* p_mti = (ci.va).c_cf ;
109
110 assert( ci.check_dzpuis(4) ) ; // dzpuis must be equal to 4
111
112 assert(p_mti->get_etat() == ETATQCQ) ;
113
114 resu->set_etat_qcq() ; // Allocates the memory for the array of double
115
116 // Loop of the various domains
117 // ---------------------------
118 for (int l=0 ; l<nz ; l++) {
119
120 const Tbl* p_tbi = p_mti->t[l] ;
121
122 if ( p_tbi->get_etat() == ETATZERO ) {
123 resu->t[l] = 0 ;
124 }
125 else { // Case where the computation must be done
126
127 assert( p_tbi->get_etat() == ETATQCQ ) ;
128
129 int nt = mg->get_nt(l) ;
130 int nr = mg->get_nr(l) ;
131
132 int base = (p_mti->base).get_b(l) ;
133 int base_r = base & MSQ_R ;
134 int base_t = base & MSQ_T ;
135 int base_p = base & MSQ_P ;
136
137 // ----------------------------------
138 // Integration on theta -> array in r
139 // ----------------------------------
140
141 double* s_tr = new double[nr] ; // Partial integral on theta
142 double* x_spec = p_tbi->t ; // Pointer on the spectral coefficients
143
144 switch (base_t) {
145
146 case T_COS_P: case T_COSSIN_CP: {
147 if (nt > nt_cp_pre) { // Initialization of factors for summation
148 nt_cp_pre = nt ;
149 cx_tcp
150 = static_cast<double*>(realloc(cx_tcp, nt*sizeof(double))) ;
151 for (int j=0 ; j<nt ; j++) {
152 cx_tcp[j] = 2./(1. - 4.*j*j) ; // Factor 2 symmetry
153 }
154 }
155
156 // Summation :
157 for (int i=0 ; i<nr ; i++) s_tr[i] = 0 ;
158 for (int j=0 ; j<nt ; j++) {
159 for (int i=0 ; i<nr ; i++) {
160 s_tr[i] += cx_tcp[j] * x_spec[i] ;
161 }
162 x_spec += nr ; // Next theta
163 }
164 break ;
165 }
166 case T_COSSIN_C: case T_COS: {
167 // Summation :
168 for (int i=0 ; i<nr ; i++) s_tr[i] = 0 ;
169 for (int j=0 ; j<nt ; j++) {
170 if ((j%2)==0)
171 for (int i=0 ; i<nr ; i++) {
172 s_tr[i] += (2. / (1.-j*j)) * x_spec[i] ;
173 }
174 x_spec += nr ; // Next theta
175 }
176 break ;
177 }
178 default: {
179 cout << "Map_af::integrale: unknown theta basis ! " << endl ;
180 abort () ;
181 break ;
182 }
183
184 } // End of the various cases on base_t
185
186 // ----------------
187 // Integration on r
188 // ----------------
189
190 double som = 0 ;
191 double som_x2 = 0;
192 double som_x = 0 ;
193 double som_c = 0 ;
194
195 switch(base_r) {
196 case R_CHEBP: case R_CHEBPIM_P: case R_CHEBPI_P :{
197 assert(beta[l] == 0) ;
198 if (nr > nr_p_pre) { // Initialization of factors for summation
199 nr_p_pre = nr ;
200 cx_rrp = static_cast<double*>(realloc(cx_rrp, nr*sizeof(double))) ;
201 for (int i=0 ; i<nr ; i++) {
202 cx_rrp[i] = (3. - 4.*i*i) /
203 (9. - 40. * i*i + 16. * i*i*i*i) ;
204 }
205 }
206
207 for (int i=0 ; i<nr ; i++) {
208 som += cx_rrp[i] * s_tr[i] ;
209 }
210 double rmax = alpha[l] ;
211 som *= rmax*rmax*rmax ;
212 break ;
213 }
214
215 case R_CHEB: {
216 if (nr > nr_f_pre) { // Initialization of factors for summation
217 nr_f_pre = nr ;
218 cx_rf_x2 = static_cast<double*>(realloc(cx_rf_x2, nr*sizeof(double))) ;
219 cx_rf_x = static_cast<double*>(realloc(cx_rf_x, nr*sizeof(double))) ;
220 cx_rf = static_cast<double*>(realloc(cx_rf, nr*sizeof(double))) ;
221 for (int i=0 ; i<nr ; i +=2 ) {
222 cx_rf_x2[i] = 2.*(3. - i*i)/(9. - 10. * i*i + i*i*i*i) ;
223 cx_rf_x[i] = 0 ;
224 cx_rf[i] = 2./(1. - i*i) ;
225 }
226 for (int i=1 ; i<nr ; i +=2 ) {
227 cx_rf_x2[i] = 0 ;
228 cx_rf_x[i] = 2./(4. - i*i) ;
229 cx_rf[i] = 0 ;
230 }
231 }
232
233 for (int i=0 ; i<nr ; i +=2 ) {
234 som_x2 += cx_rf_x2[i] * s_tr[i] ;
235 }
236 for (int i=1 ; i<nr ; i +=2 ) {
237 som_x += cx_rf_x[i] * s_tr[i] ;
238 }
239 for (int i=0 ; i<nr ; i +=2 ) {
240 som_c += cx_rf[i] * s_tr[i] ;
241 }
242 double a = alpha[l] ;
243 double b = beta[l] ;
244 som = a*a*a * som_x2 + 2.*a*a*b * som_x + a*b*b * som_c ;
245 break ;
246 }
247
248 case R_JACO02: {
249 if (nr > nr_f_pre) { // Initialization of factors for summation
250 nr_f_pre = nr ;
251 cx_rf_x2 = static_cast<double*>(realloc(cx_rf_x2, nr*sizeof(double))) ;
252 cx_rf_x = static_cast<double*>(realloc(cx_rf_x, nr*sizeof(double))) ;
253 cx_rf = static_cast<double*>(realloc(cx_rf, nr*sizeof(double))) ;
254 double signe = 1 ;
255 for (int i=0 ; i<nr ; i +=1 ) {
256 cx_rf_x2[i] = 0 ;
257 cx_rf_x[i] = 2*signe/double(i+1)/double(i+2);
258 cx_rf[i] = 2*signe ;
259 signe = - signe ;
260 }
261 cx_rf_x2[0] = double(8)/double(3) ;
262 }
263
264 for (int i=0 ; i<nr ; i +=1 ) {
265 som_x2 += cx_rf_x2[i] * s_tr[i] ;
266 }
267 for (int i=1 ; i<nr ; i +=1 ) {
268 som_x += cx_rf_x[i] * s_tr[i] ;
269 }
270 for (int i=0 ; i<nr ; i +=1 ) {
271 som_c += cx_rf[i] * s_tr[i] ;
272 }
273 double a = alpha[l] ;
274 double b = beta[l] ;
275 assert(b == a) ;
276 som = a*a*a * som_x2 + 2.*a*a*(b-a) * som_x + a*(b-a)*(b-a) * som_c ;
277 break ;
278 }
279
280 case R_CHEBU: {
281 assert(beta[l] == - alpha[l]) ;
282 if (nr > nr_f_pre) { // Initialization of factors for summation
283 nr_f_pre = nr ;
284 cx_rf = static_cast<double*>(realloc(cx_rf, nr*sizeof(double))) ;
285 for (int i=0 ; i<nr ; i +=2 ) {
286 cx_rf[i] = 2./(1. - i*i) ;
287 }
288 for (int i=1 ; i<nr ; i +=2 ) {
289 cx_rf[i] = 0 ;
290 }
291 }
292
293 for (int i=0 ; i<nr ; i +=2 ) {
294 som_c += cx_rf[i] * s_tr[i] ;
295 }
296 som = - alpha[l] * som_c ;
297 break ;
298 }
299
300
301 default: {
302 cout << "Map_af::integrale: unknown r basis ! " << endl ;
303 abort () ;
304 break ;
305 }
306
307 } // End of the various cases on base_r
308
309 // ------------------
310 // Integration on phi
311 // ------------------
312
313 switch (base_p) {
314
315 case P_COSSIN: {
316 som *= 2. * M_PI ;
317 break ;
318 }
319
320 case P_COSSIN_P: {
321 som *= 2. * M_PI ;
322 break ;
323 }
324
325 default: {
326 cout << "Map_af::integrale: unknown phi basis ! " << endl ;
327 abort () ;
328 break ;
329 }
330
331 } // End of the various cases on base_p
332
333 // Final result for this domain:
334 // ----------------------------
335
336 resu->t[l] = som ;
337 delete [] s_tr ;
338
339 } // End of the case where the computation must be done
340
341
342 } // End of the loop onto the domains
343
344 return resu ;
345
346}
347}
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
int get_etat() const
Returns the logical state.
Definition cmp.h:899
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
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 cmp.C:715
double * beta
Array (size: mg->nzone ) of the values of in each domain.
Definition map.h:2035
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition map.h:2033
virtual Tbl * integrale(const Cmp &) const
Computes the integral over all space of a Cmp.
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition map.h:676
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
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
int get_etat() const
Returns the logical state.
Definition mtbl_cf.h:456
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
Definition mtbl_cf.h:205
Basic array class.
Definition tbl.h:161
int get_etat() const
Gives the logical state.
Definition tbl.h:394
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:372
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double * t
The array of double.
Definition tbl.h:173
#define P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
#define MSQ_R
Extraction de l'info sur R.
#define T_COS_P
dev. cos seulement, harmoniques paires
#define P_COSSIN
dev. standart
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
#define MSQ_T
Extraction de l'info sur Theta.
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
#define R_CHEB
base de Chebychev ordinaire (fin)
#define T_COS
dev. cos seulement
#define R_CHEBP
base de Cheb. paire (rare) seulement
#define MSQ_P
Extraction de l'info sur Phi.
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Lorene prototypes.
Definition app_hor.h:64