LORENE
scalar_raccord_externe.C
1 /*
2  * Copyright (c) 2003 Eric Gourgoulhon & Jerome Novak
3  *
4  * Copyright (c) 2001 Philippe Grandclement (for preceding Cmp version)
5  *
6  *
7  * This file is part of LORENE.
8  *
9  * LORENE is free software; you can redistribute it and/or modify
10  * it under the terms of the GNU General Public License as published by
11  * the Free Software Foundation; either version 2 of the License, or
12  * (at your option) any later version.
13  *
14  * LORENE is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the GNU General Public License
20  * along with LORENE; if not, write to the Free Software
21  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22  *
23  */
24 
25 
26 char scalar_raccord_externe_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_raccord_externe.C,v 1.4 2014/10/13 08:53:47 j_novak Exp $" ;
27 
28 /*
29  * $Id: scalar_raccord_externe.C,v 1.4 2014/10/13 08:53:47 j_novak Exp $
30  * $Log: scalar_raccord_externe.C,v $
31  * Revision 1.4 2014/10/13 08:53:47 j_novak
32  * Lorene classes and functions now belong to the namespace Lorene.
33  *
34  * Revision 1.3 2014/10/06 15:16:16 j_novak
35  * Modified #include directives to use c++ syntax.
36  *
37  * Revision 1.2 2003/09/25 09:22:33 j_novak
38  * Added a #include<math.h>
39  *
40  * Revision 1.1 2003/09/25 08:58:10 e_gourgoulhon
41  * First version.
42  *
43  *
44  * $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_raccord_externe.C,v 1.4 2014/10/13 08:53:47 j_novak Exp $
45  *
46  */
47 
48 
49 
50 //standard
51 #include <cstdlib>
52 #include <cmath>
53 
54 // LORENE
55 #include "matrice.h"
56 #include "tensor.h"
57 #include "proto.h"
58 
59 namespace Lorene {
60 int cnp (int n, int p) ;
61 
62 // Fait le raccord dans la zec ...
63 // Suppose (pour le moment, le meme nbre de points sur les angles ...)
64 // et que la zone precedente est une coquille
65 
66 void Scalar::raccord_externe(int power, int nbre, int lmax) {
67 
68  va.coef() ;
69  va.ylm() ;
70 
71  Base_val base_devel (va.base) ;
72  int base_r, m_quant, l_quant ;
73 
74  // Confort :
75  int zone = mp->get_mg()->get_nzone()-2 ;
76  int nt = mp->get_mg()->get_nt(zone) ;
77  int np = mp->get_mg()->get_np(zone) ;
78  int nr = mp->get_mg()->get_nr(zone) ;
79 
80  // Le mapping doit etre affine :
81  const Map_af* map = dynamic_cast<const Map_af*>(mp) ;
82  if (map == 0x0) {
83  cout << "Le mapping doit etre affine" << endl ;
84  abort() ;
85  }
86 
87  // Mappinhg en r
88  double alpha = map->get_alpha()[zone] ;
89  double beta = map->get_beta()[zone] ;
90 
91  // Mapping en 1/r
92  double new_alpha = -alpha/(beta*beta-alpha*alpha) ;
93  double new_beta = beta/(beta*beta-alpha*alpha) ;
94 
95  // Mapping dans la zec :
96  double alpha_zec = map->get_alpha()[zone+1] ;
97 
98  // Maintenant on construit les matrices de passage :
99  // Celle de ksi a T
100  Matrice tksi (nbre, nbre) ;
101  tksi.set_etat_qcq() ;
102 
103  // Premier polynome
104  tksi.set(0, 0) = sqrt(double(2)) ;
105  for (int i=1 ; i<nbre ; i++)
106  tksi.set(0, i) = 0 ;
107 
108  //Second polynome
109  tksi.set(1, 0) = 0 ;
110  tksi.set(1, 1) = sqrt(double(2)) ;
111  for (int i=2 ; i<nbre ; i++)
112  tksi.set(1, i) = 0 ;
113 
114  // On recurre :
115  for (int lig=2 ; lig<nbre ; lig++) {
116  tksi.set(lig, 0) = -tksi(lig-2, 0) ;
117  for (int col=1 ; col<nbre ; col++)
118  tksi.set(lig, col) = 2*tksi(lig-1, col-1)-tksi(lig-2, col) ;
119  }
120 
121  // Celle de u/new_alpha a ksi :
122  Matrice ksiu (nbre, nbre) ;
123  ksiu.set_etat_qcq() ;
124 
125  for (int lig=0 ; lig<nbre ; lig++) {
126  for (int col=0 ; col<=lig ; col++)
127  ksiu.set(lig, col) = cnp(lig, col)*
128  pow(-new_beta/new_alpha, lig-col) ;
129  for (int col = lig+1 ; col<nbre ; col++)
130  ksiu.set(lig, col) = 0 ;
131  }
132 
133  // La matrice totale :
134  Matrice tu (nbre, nbre) ;
135  tu.set_etat_qcq() ;
136  double somme ;
137  for (int lig=0 ; lig<nbre ; lig++)
138  for (int col=0 ; col<nbre ; col++) {
139  somme = 0 ;
140  for (int m=0 ; m<nbre ; m++)
141  somme += tksi(lig, m)*ksiu(m, col) ;
142  tu.set(lig, col) = somme ;
143  }
144 
145  // On calcul les coefficients de u^n dans la zec
146  Tbl coef_u (nbre+lmax, nr) ;
147  coef_u.set_etat_qcq() ;
148  int* dege = new int [3] ;
149  dege[0] = 1 ; dege[1] = 1 ; dege[2] = nr ;
150  double* ti = new double [nr] ;
151 
152  for (int puiss=0 ; puiss<nbre+lmax ; puiss++) {
153  for (int i=0 ; i<nr ; i++)
154  ti[i] = pow(-cos(M_PI*i/(nr-1))-1, puiss) ;
155  cfrcheb (dege, dege, ti, dege, ti) ;
156  for (int i=0 ; i<nr ; i++)
157  coef_u.set(puiss, i) = ti[i] ;
158  }
159 
160  // Avant d entrer dans la boucle :
161  dege[2] = nbre ;
162  double *coloc = new double[nbre] ;
163  double *auxi = new double [1] ;
164 
165  Tbl coef_zec (np+2, nt, nr) ;
166  coef_zec.annule_hard() ;
167 
168  // Boucle sur les harmoniques :
169 
170  for (int k=0 ; k<np+2 ; k++)
171  for (int j=0 ; j<nt ; j++)
172  if (nullite_plm (j, nt, k, np, base_devel)==1) {
173  donne_lm (zone+2, zone+1, j, k, base_devel, m_quant,
174  l_quant, base_r) ;
175  if (l_quant <= lmax) {
176 
177  // On bosse :
178  // On recupere les valeus aux points de colocation en 1/r :
179  double ksi, air ;
180  for (int i=0 ; i<nbre ; i++) {
181  ksi = -cos(M_PI*i/(nbre-1)) ;
182  air = 1./(new_alpha*ksi+new_beta) ;
183  ksi = (air-beta)/alpha ;
184  for (int m=0 ; m<nr ; m++)
185  ti[m] = (*va.c_cf)(zone, k, j, m) ;
186  som_r_cheb (ti, nr, 1, 1, ksi, auxi) ;
187  coloc[i] = auxi[0]/
188  pow (-new_alpha*cos(M_PI*i/(nbre-1))+new_beta, power+l_quant);
189  }
190 
191  cfrcheb (dege, dege, coloc, dege, coloc) ;
192 
193  Tbl expansion (nbre) ;
194  expansion.set_etat_qcq() ;
195  for (int i=0 ; i<nbre ; i++) {
196  somme = 0 ;
197  for (int m=0 ; m<nbre ; m++)
198  somme += coloc[m]*tu(m, i) ;
199  expansion.set(i) = somme ;
200  }
201 
202  for (int i=0 ; i<nr ; i++) {
203  somme = 0 ;
204  for (int m=0 ; m<nbre ; m++)
205  somme += coef_u(m+l_quant, i)*expansion(m)*
206  pow(alpha_zec, m+l_quant)/
207  pow(new_alpha, m) ;
208  coef_zec.set(k, j, i) = somme ;
209  }
210  }
211  }
212 
213  va.set_etat_cf_qcq() ;
214  va.c_cf->set_etat_qcq() ;
215  va.c_cf->t[zone+1]->set_etat_qcq() ;
216 
217  for (int k=0 ; k<np+2 ; k++)
218  for (int j=0 ; j<nt ; j++)
219  for (int i=0 ; i<nr ; i++)
220  va.c_cf->set(zone+1, k, j, i) = coef_zec(k, j, i) ;
221 
222  set_dzpuis(power) ;
223  va.ylm_i() ;
224 
225  delete[] auxi ;
226  delete [] dege ;
227  delete [] ti ;
228  delete [] coloc ;
229 }
230 }
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Lorene::Base_val
Bases of the spectral expansions.
Definition: base_val.h:322
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Lorene::Valeur::set_etat_cf_qcq
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
Lorene::Matrice::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:175
Lorene::Scalar::va
Valeur va
The numerical value of the Scalar
Definition: scalar.h:405
Lorene::Mtbl_cf::t
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
Definition: mtbl_cf.h:205
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:372
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Scalar::raccord_externe
void raccord_externe(int puis, int nbre, int lmax)
Matching of the external domain with the outermost shell.
Definition: scalar_raccord_externe.C:66
Lorene::Mtbl_cf::set
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:294
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Lorene::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:295
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Lorene::Scalar::sqrt
friend Scalar sqrt(const Scalar &)
Square root.
Definition: scalar_math.C:263
Lorene::Map_af
Affine radial mapping.
Definition: map.h:2027
Lorene::Scalar::set_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:808
Lorene::Matrice
Matrix handling.
Definition: matrice.h:152
Lorene::Mtbl_cf::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl_cf.C:300
Lorene::Map_af::get_beta
const double * get_beta() const
Returns the pointer on the array beta.
Definition: map_af.C:481
Lorene::Valeur::base
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
Lorene::Matrice::set
double & set(int j, int i)
Read/write of a particuliar element.
Definition: matrice.h:277
Lorene::Scalar::pow
friend Scalar pow(const Scalar &, int)
Power .
Definition: scalar_math.C:451
Lorene::Map_af::get_alpha
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:477
Lorene::Scalar::cos
friend Scalar cos(const Scalar &)
Cosine.
Definition: scalar_math.C:104
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131