LORENE
map_af_integ_surf.C
1 /*
2  * Copyright (c) 2000-2001 Philippe Grandclement
3  *
4  * This file is part of LORENE.
5  *
6  * LORENE is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * LORENE is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with LORENE; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 
23 char map_af_integ_surf_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_af_integ_surf.C,v 1.7 2014/10/13 08:53:02 j_novak Exp $" ;
24 
25 /*
26  * $Id: map_af_integ_surf.C,v 1.7 2014/10/13 08:53:02 j_novak Exp $
27  * $Log: map_af_integ_surf.C,v $
28  * Revision 1.7 2014/10/13 08:53:02 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.6 2014/10/06 15:13:12 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.5 2009/10/08 16:20:47 j_novak
35  * Addition of new bases T_COS and T_SIN.
36  *
37  * Revision 1.4 2007/10/05 15:56:19 j_novak
38  * Addition of new bases for the non-symmetric case in theta.
39  *
40  * Revision 1.3 2004/03/10 12:43:06 jl_jaramillo
41  * Treatment of case ETATUN in surface integrals for Scalar's.
42  *
43  * Revision 1.2 2004/01/29 08:50:03 p_grandclement
44  * Modification of Map::operator==(const Map&) and addition of the surface
45  * integrales using Scalar.
46  *
47  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
48  * LORENE
49  *
50  * Revision 1.7 2001/02/19 11:40:27 phil
51  * correction indices
52  *
53  * Revision 1.6 2001/02/12 14:14:29 phil
54  * *** empty log message ***
55  *
56  * Revision 1.5 2001/02/12 14:00:52 phil
57  * on prends tous les coefficients now
58  *
59  * Revision 1.4 2001/02/12 12:35:34 phil
60  * gestion des bases angulaires plus proprement
61  *
62  * Revision 1.3 2001/01/02 10:52:27 phil
63  * ajout calcul a l'infini
64  *
65  * Revision 1.2 2000/09/19 13:54:14 phil
66  * *** empty log message ***
67  *
68  * Revision 1.1 2000/09/19 13:09:32 phil
69  * Initial revision
70  *
71  *
72  * $Header: /cvsroot/Lorene/C++/Source/Map/map_af_integ_surf.C,v 1.7 2014/10/13 08:53:02 j_novak Exp $
73  *
74  */
75 
76 
77 #include <cstdlib>
78 #include <cmath>
79 
80 #include "map.h"
81 #include "cmp.h"
82 #include "proto.h"
83 #include "scalar.h"
84 
85  //=============
86  // Cmp version
87  //===============
88 
89 namespace Lorene {
90 double Map_af::integrale_surface (const Cmp& ci, double rayon) const {
91 
92  assert (ci.get_etat() != ETATNONDEF) ;
93  assert (rayon > 0) ;
94  if (ci.get_etat() == ETATZERO)
95  return 0 ;
96 
97  assert (ci.get_etat() == ETATQCQ) ;
98 
99  int l ;
100  double xi ;
101  val_lx (rayon, 0, 0, l, xi) ;
102 
103  if (l == get_mg()->get_nzone()-1) {
104  ci.check_dzpuis(0) ;
105  }
106 
107  ci.va.coef() ;
108  int nr = get_mg()->get_nr(l) ;
109  int nt = get_mg()->get_nt(l) ;
110 
111  int base_r = ci.va.base.get_base_r(l) ;
112  int base_t = ci.va.base.get_base_t(l) ;
113  int base_p = ci.va.base.get_base_p(l) ;
114 
115  double result = 0 ;
116  double* coef = new double [nr] ;
117  double* auxi = new double[1] ;
118 
119  bool odd_theta = false ;
120  double c_cos = 2 ;
121  switch (base_t) {
122  case T_COS_P : case T_COSSIN_CP :
123  break ;
124  case T_COS_I : case T_COSSIN_CI :
125  odd_theta = true ;
126  break ;
127  case T_COS : case T_COSSIN_C :
128  c_cos = 1. ;
129  break ;
130  default :
131  cout << "base_t cas non prevu dans Map_af::integrale_surface" << endl ;
132  abort() ;
133  break ;
134  }
135 
136  if (!odd_theta) {
137  for (int j=0 ; j<nt ; j++) {
138  for (int i=0 ; i<nr ; i++)
139  coef[i] = (*ci.va.c_cf)(l, 0, j, i) ;
140 
141  switch (base_r) {
142 
143  case R_CHEB :
144  som_r_cheb (coef, nr, 1, 1, xi, auxi) ;
145  break ;
146  case R_CHEBP :
147  som_r_chebp (coef, nr, 1, 1, xi, auxi) ;
148  break ;
149  case R_CHEBI :
150  som_r_chebi (coef, nr, 1, 1, xi, auxi) ;
151  break ;
152  case R_CHEBU :
153  som_r_chebu (coef, nr, 1, 1, xi, auxi) ;
154  break ;
155  case R_CHEBPI_P :
156  som_r_chebpi_p (coef, nr, 1, 1, xi, auxi) ;
157  break ;
158  case R_CHEBPI_I :
159  som_r_chebpi_i (coef, nr, 1, 1, xi, auxi) ;
160  break ;
161  case R_CHEBPIM_P :
162  som_r_chebpim_p (coef, nr, 1, 1, xi, auxi) ;
163  break ;
164  case R_CHEBPIM_I :
165  som_r_chebpim_i (coef, nr, 1, 1, xi, auxi) ;
166  break ;
167  default :
168  som_r_pas_prevu (coef, nr, 1, 1, xi, auxi) ;
169  break ;
170  }
171  result += 2 * (*auxi)/(1-c_cos*c_cos*j*j) ;
172  if (c_cos == 1.) j++ ;
173  }
174  }
175  delete [] auxi ;
176  delete [] coef ;
177 
178  switch (base_p) {
179  case P_COSSIN :
180  result *= 2*rayon*rayon*M_PI ;
181  break ;
182  case P_COSSIN_P :
183  result *= 2*rayon*rayon*M_PI ;
184  break ;
185  default :
186  cout << "base_p cas non prevu dans Map_af::integrale_surface" << endl ;
187  abort() ;
188  break ;
189  }
190 
191  return result ;
192 }
193 
194 
195 // Integrale a l'infini
196 double Map_af::integrale_surface_infini (const Cmp& ci) const {
197 
198  assert (ci.get_etat() != ETATNONDEF) ;
199  assert (ci.check_dzpuis(2));
200 
201  if (ci.get_etat() == ETATZERO)
202  return 0 ;
203 
204  assert (ci.get_etat() == ETATQCQ) ;
205 
206  int nz = ci.get_mp()->get_mg()->get_nzone() ;
207 
208  ci.va.coef() ;
209  int nr = get_mg()->get_nr(nz-1) ;
210  int nt = get_mg()->get_nt(nz-1) ;
211 
212  int base_r = ci.va.base.get_base_r(nz-1) ;
213  int base_t = ci.va.base.get_base_t(nz-1) ;
214  int base_p = ci.va.base.get_base_p(nz-1) ;
215 
216  double result = 0 ;
217  double* coef = new double [nr] ;
218  double* auxi = new double[1] ;
219 
220  bool odd_theta = false ;
221  double c_cos = 2. ;
222  switch (base_t) {
223  case T_COS_P : case T_COSSIN_CP :
224  break ;
225  case T_COS_I : case T_COSSIN_CI :
226  odd_theta = true ;
227  break ;
228  case T_COS : case T_COSSIN_C :
229  c_cos = 1. ;
230  break ;
231  default :
232  cout << "base_t cas non prevu dans Map_af::integrale_surface" << endl ;
233  abort() ;
234  break ;
235  }
236 
237  if (!odd_theta) {
238  for (int j=0 ; j<nt ; j++) {
239  for (int i=0 ; i<nr ; i++)
240  coef[i] = (*ci.va.c_cf)(nz-1, 0, j, i) ;
241 
242  switch (base_r) {
243  case R_CHEBU :
244  som_r_chebu (coef, nr, 1, 1, 1, auxi) ;
245  break ;
246  default :
247  som_r_pas_prevu (coef, nr, 1, 1, 1, auxi) ;
248  break ;
249  }
250  result += 2 * (*auxi)/(1-c_cos*c_cos*j*j) ;
251  if (c_cos == 1.) j++ ;
252  }
253  }
254  delete [] auxi ;
255  delete [] coef ;
256 
257  switch (base_p) {
258  case P_COSSIN :
259  result *= 2*M_PI ;
260  break ;
261  case P_COSSIN_P :
262  result *= 2*M_PI ;
263  break ;
264  default :
265  cout << "base_p cas non prevu dans Map_af::integrale_surface_infini" << endl ;
266  abort() ;
267  break ;
268  }
269 
270  return result ;
271 }
272 
273  //=============
274  // Scalar version
275  //===============
276 
277 double Map_af::integrale_surface (const Scalar& ci, double rayon) const {
278 
279  assert (ci.get_etat() != ETATNONDEF) ;
280  assert (rayon > 0) ;
281  if (ci.get_etat() == ETATZERO)
282  return 0 ;
283 
284  assert ( (ci.get_etat() == ETATQCQ) || (ci.get_etat() == ETATUN) ) ;
285 
286  int l ;
287  double xi ;
288  val_lx (rayon, 0, 0, l, xi) ;
289 
290  if (l == get_mg()->get_nzone()-1) {
291  ci.check_dzpuis(0) ;
292  }
293 
294  ci.get_spectral_va().coef() ;
295  int nr = get_mg()->get_nr(l) ;
296  int nt = get_mg()->get_nt(l) ;
297 
298  int base_r = ci.get_spectral_va().base.get_base_r(l) ;
299  int base_t = ci.get_spectral_va().base.get_base_t(l) ;
300  int base_p = ci.get_spectral_va().base.get_base_p(l) ;
301 
302  double result = 0 ;
303  double* coef = new double [nr] ;
304  double* auxi = new double[1] ;
305 
306  bool odd_theta = false ;
307  double c_cos = 2. ;
308 
309  switch (base_t) {
310  case T_COS_P : case T_COSSIN_CP :
311  break ;
312  case T_COS_I: case T_COSSIN_CI :
313  odd_theta = true ;
314  break ;
315  case T_COS : case T_COSSIN_C :
316  c_cos = 1. ;
317  break ;
318  default :
319  cout << "base_t cas non prevu dans Map_af::integrale_surface" << endl ;
320  abort() ;
321  break ;
322  }
323 
324  if (!odd_theta) {
325  for (int j=0 ; j<nt ; j++) {
326  for (int i=0 ; i<nr ; i++)
327  coef[i] = (*ci.get_spectral_va().c_cf)(l, 0, j, i) ;
328 
329  switch (base_r) {
330 
331  case R_CHEB :
332  som_r_cheb (coef, nr, 1, 1, xi, auxi) ;
333  break ;
334  case R_CHEBP :
335  som_r_chebp (coef, nr, 1, 1, xi, auxi) ;
336  break ;
337  case R_CHEBI :
338  som_r_chebi (coef, nr, 1, 1, xi, auxi) ;
339  break ;
340  case R_CHEBU :
341  som_r_chebu (coef, nr, 1, 1, xi, auxi) ;
342  break ;
343  case R_CHEBPI_P :
344  som_r_chebpi_p (coef, nr, 1, 1, xi, auxi) ;
345  break ;
346  case R_CHEBPI_I :
347  som_r_chebpi_i (coef, nr, 1, 1, xi, auxi) ;
348  break ;
349  case R_CHEBPIM_P :
350  som_r_chebpim_p (coef, nr, 1, 1, xi, auxi) ;
351  break ;
352  case R_CHEBPIM_I :
353  som_r_chebpim_i (coef, nr, 1, 1, xi, auxi) ;
354  break ;
355  default :
356  som_r_pas_prevu (coef, nr, 1, 1, xi, auxi) ;
357  break ;
358  }
359  result += 2 * (*auxi)/(1-c_cos*c_cos*j*j) ;
360  if (c_cos == 1.) j++ ;
361  }
362  }
363 
364  delete [] auxi ;
365  delete [] coef ;
366 
367  switch (base_p) {
368  case P_COSSIN :
369  result *= 2*rayon*rayon*M_PI ;
370  break ;
371  case P_COSSIN_P :
372  result *= 2*rayon*rayon*M_PI ;
373  break ;
374  default :
375  cout << "base_p cas non prevu dans Map_af::integrale_surface" << endl ;
376  abort() ;
377  break ;
378  }
379 
380  return result ;
381 }
382 
383 
384 // Integrale a l'infini
385 double Map_af::integrale_surface_infini (const Scalar& ci) const {
386 
387  assert (ci.get_etat() != ETATNONDEF) ;
388  assert (ci.check_dzpuis(2));
389 
390  if (ci.get_etat() == ETATZERO)
391  return 0 ;
392 
393  assert ( (ci.get_etat() == ETATQCQ) || (ci.get_etat() == ETATUN) ) ;
394 
395  int nz = ci.get_mp().get_mg()->get_nzone() ;
396 
397  ci.get_spectral_va().coef() ;
398  int nr = get_mg()->get_nr(nz-1) ;
399  int nt = get_mg()->get_nt(nz-1) ;
400 
401  int base_r = ci.get_spectral_va().base.get_base_r(nz-1) ;
402  int base_t = ci.get_spectral_va().base.get_base_t(nz-1) ;
403  int base_p = ci.get_spectral_va().base.get_base_p(nz-1) ;
404 
405  double result = 0 ;
406  double* coef = new double [nr] ;
407  double* auxi = new double[1] ;
408 
409  bool odd_theta = false ;
410  double c_cos = 2. ;
411 
412  switch (base_t) {
413  case T_COS_P : case T_COSSIN_CP :
414  break ;
415  case T_COS_I: case T_COSSIN_CI :
416  odd_theta = true ;
417  break ;
418  case T_COS : case T_COSSIN_C :
419  c_cos = 1. ;
420  break ;
421  default :
422  cout << "base_t cas non prevu dans Map_af::integrale_surface" << endl ;
423  abort() ;
424  break ;
425  }
426 
427  if (!odd_theta) {
428  for (int j=0 ; j<nt ; j++) {
429  for (int i=0 ; i<nr ; i++)
430  coef[i] = (*ci.get_spectral_va().c_cf)(nz-1, 0, j, i) ;
431 
432  switch (base_r) {
433  case R_CHEBU :
434  som_r_chebu (coef, nr, 1, 1, 1, auxi) ;
435  break ;
436  default :
437  som_r_pas_prevu (coef, nr, 1, 1, 1, auxi) ;
438  break ;
439  }
440  result += 2 * (*auxi)/(1-c_cos*c_cos*j*j) ;
441  if (c_cos == 1.) j++ ;
442  }
443  }
444 
445  delete [] auxi ;
446  delete [] coef ;
447 
448  switch (base_p) {
449  case P_COSSIN :
450  result *= 2*M_PI ;
451  break ;
452  case P_COSSIN_P :
453  result *= 2*M_PI ;
454  break ;
455  default :
456  cout << "base_p cas non prevu dans Map_af::integrale_surface_infini" << endl ;
457  abort() ;
458  break ;
459  }
460 
461  return result ;
462 }
463 }
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
#define P_COSSIN
dev. standart
Definition: type_parite.h:245
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
double integrale_surface_infini(const Cmp &ci) const
Performs the surface integration of ci at infinity.
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition: type_parite.h:174
Lorene prototypes.
Definition: app_hor.h:64
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
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
#define T_COS
dev. cos seulement
Definition: type_parite.h:196
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
int get_base_r(int l) const
Returns the expansion basis for r ( ) functions in the domain of index l (e.g.
Definition: base_val.h:400
#define T_COS_I
dev. cos seulement, harmoniques impaires
Definition: type_parite.h:204
#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
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition: type_parite.h:200
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition: type_parite.h:192
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
int get_base_p(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
Definition: base_val.h:422
virtual void val_lx(double rr, double theta, double pphi, int &l, double &xi) const
Computes the domain index l and the value of corresponding to a point given by its physical coordina...
#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
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition: type_parite.h:172
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
Definition: type_parite.h:212
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
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
#define P_COSSIN_P
dev. sur Phi = 2*phi, freq. paires
Definition: type_parite.h:247
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
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
Definition: type_parite.h:208
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:861
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:601
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166