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 
28 char 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 
80 namespace Lorene {
81 Tbl* 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.
Definition: map_af_integ.C:81
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
Definition: type_parite.h:247
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
#define MSQ_R
Extraction de l'info sur R.
Definition: type_parite.h:152
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition: type_parite.h:200
#define P_COSSIN
dev. standart
Definition: type_parite.h:245
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
Definition: type_parite.h:176
#define MSQ_T
Extraction de l'info sur Theta.
Definition: type_parite.h:154
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
Definition: type_parite.h:208
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
#define T_COS
dev. cos seulement
Definition: type_parite.h:196
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
#define MSQ_P
Extraction de l'info sur Phi.
Definition: type_parite.h:156
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition: type_parite.h:192
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition: type_parite.h:172
Lorene prototypes.
Definition: app_hor.h:64