LORENE
cmp_raccord_externe.C
1 /*
2  * Copyright (c) 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 cmp_raccord_externe_C[] = "$Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_raccord_externe.C,v 1.4 2014/10/13 08:52:48 j_novak Exp $" ;
24 
25 /*
26  * $Id: cmp_raccord_externe.C,v 1.4 2014/10/13 08:52:48 j_novak Exp $
27  * $Log: cmp_raccord_externe.C,v $
28  * Revision 1.4 2014/10/13 08:52:48 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.3 2014/10/06 15:13:04 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.2 2003/10/03 15:58:45 j_novak
35  * Cleaning of some headers
36  *
37  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
38  * LORENE
39  *
40  * Revision 2.2 2001/10/10 13:53:27 eric
41  * Modif Joachim: sqrt(2) --> sqrt(double(2))
42  *
43  * Revision 2.1 2001/04/02 12:16:39 phil
44  * *** empty log message ***
45  *
46  * Revision 2.0 2001/03/30 13:37:32 phil
47  * *** empty log message ***
48  *
49  *
50  * $Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_raccord_externe.C,v 1.4 2014/10/13 08:52:48 j_novak Exp $
51  *
52  */
53 
54 
55 
56 //standard
57 #include <cstdlib>
58 #include <cmath>
59 
60 // LORENE
61 #include "matrice.h"
62 #include "cmp.h"
63 #include "proto.h"
64 
65 
66 // Calcul des Cnp
67 namespace Lorene {
68 int cnp (int n, int p) {
69 
70  assert (p<=n) ;
71 
72  if ((p==0) || (p==n))
73  return 1 ;
74  else {
75  int fact_un = 1 ;
76  for (int conte=n ; conte >n-p ; conte --)
77  fact_un *= conte ;
78 
79  int fact_deux = 1 ;
80  for (int conte = 1 ; conte<p+1 ; conte++)
81  fact_deux *= conte ;
82 
83  return int(fact_un/fact_deux) ;
84  }
85 }
86 
87 // Fait le raccord dans la zec ...
88 // Suppose (pour le moment, le meme nbre de points sur les angles ...)
89 // et que la zone precedente est une coquille
90 
91 void Cmp::raccord_externe (int power, int nbre, int lmax) {
92 
93  va.coef() ;
94  va.ylm() ;
95 
96  Base_val base_devel (va.base) ;
97  int base_r, m_quant, l_quant ;
98 
99  // Confort :
100  int zone = mp->get_mg()->get_nzone()-2 ;
101  int nt = mp->get_mg()->get_nt(zone) ;
102  int np = mp->get_mg()->get_np(zone) ;
103  int nr = mp->get_mg()->get_nr(zone) ;
104 
105  // Le mapping doit etre affine :
106  const Map_af* map = dynamic_cast<const Map_af*>(mp) ;
107  if (map == 0x0) {
108  cout << "Le mapping doit etre affine" << endl ;
109  abort() ;
110  }
111 
112  // Mappinhg en r
113  double alpha = map->get_alpha()[zone] ;
114  double beta = map->get_beta()[zone] ;
115 
116  // Mapping en 1/r
117  double new_alpha = -alpha/(beta*beta-alpha*alpha) ;
118  double new_beta = beta/(beta*beta-alpha*alpha) ;
119 
120  // Mapping dans la zec :
121  double alpha_zec = map->get_alpha()[zone+1] ;
122 
123  // Maintenant on construit les matrices de passage :
124  // Celle de ksi a T
125  Matrice tksi (nbre, nbre) ;
126  tksi.set_etat_qcq() ;
127 
128  // Premier polynome
129  tksi.set(0, 0) = sqrt(double(2)) ;
130  for (int i=1 ; i<nbre ; i++)
131  tksi.set(0, i) = 0 ;
132 
133  //Second polynome
134  tksi.set(1, 0) = 0 ;
135  tksi.set(1, 1) = sqrt(double(2)) ;
136  for (int i=2 ; i<nbre ; i++)
137  tksi.set(1, i) = 0 ;
138 
139  // On recurre :
140  for (int lig=2 ; lig<nbre ; lig++) {
141  tksi.set(lig, 0) = -tksi(lig-2, 0) ;
142  for (int col=1 ; col<nbre ; col++)
143  tksi.set(lig, col) = 2*tksi(lig-1, col-1)-tksi(lig-2, col) ;
144  }
145 
146  // Celle de u/new_alpha a ksi :
147  Matrice ksiu (nbre, nbre) ;
148  ksiu.set_etat_qcq() ;
149 
150  for (int lig=0 ; lig<nbre ; lig++) {
151  for (int col=0 ; col<=lig ; col++)
152  ksiu.set(lig, col) = cnp(lig, col)*
153  pow(-new_beta/new_alpha, lig-col) ;
154  for (int col = lig+1 ; col<nbre ; col++)
155  ksiu.set(lig, col) = 0 ;
156  }
157 
158  // La matrice totale :
159  Matrice tu (nbre, nbre) ;
160  tu.set_etat_qcq() ;
161  double somme ;
162  for (int lig=0 ; lig<nbre ; lig++)
163  for (int col=0 ; col<nbre ; col++) {
164  somme = 0 ;
165  for (int m=0 ; m<nbre ; m++)
166  somme += tksi(lig, m)*ksiu(m, col) ;
167  tu.set(lig, col) = somme ;
168  }
169 
170  // On calcul les coefficients de u^n dans la zec
171  Tbl coef_u (nbre+lmax, nr) ;
172  coef_u.set_etat_qcq() ;
173  int* dege = new int [3] ;
174  dege[0] = 1 ; dege[1] = 1 ; dege[2] = nr ;
175  double* ti = new double [nr] ;
176 
177  for (int puiss=0 ; puiss<nbre+lmax ; puiss++) {
178  for (int i=0 ; i<nr ; i++)
179  ti[i] = pow(-cos(M_PI*i/(nr-1))-1, puiss) ;
180  cfrcheb (dege, dege, ti, dege, ti) ;
181  for (int i=0 ; i<nr ; i++)
182  coef_u.set(puiss, i) = ti[i] ;
183  }
184 
185  // Avant d entrer dans la boucle :
186  dege[2] = nbre ;
187  double *coloc = new double[nbre] ;
188  double *auxi = new double [1] ;
189 
190  Tbl coef_zec (np+2, nt, nr) ;
191  coef_zec.annule_hard() ;
192 
193  // Boucle sur les harmoniques :
194 
195  for (int k=0 ; k<np+2 ; k++)
196  for (int j=0 ; j<nt ; j++)
197  if (nullite_plm (j, nt, k, np, base_devel)==1) {
198  donne_lm (zone+2, zone+1, j, k, base_devel, m_quant,
199  l_quant, base_r) ;
200  if (l_quant <= lmax) {
201 
202  // On bosse :
203  // On recupere les valeus aux points de colocation en 1/r :
204  double ksi, air ;
205  for (int i=0 ; i<nbre ; i++) {
206  ksi = -cos(M_PI*i/(nbre-1)) ;
207  air = 1./(new_alpha*ksi+new_beta) ;
208  ksi = (air-beta)/alpha ;
209  for (int m=0 ; m<nr ; m++)
210  ti[m] = (*va.c_cf)(zone, k, j, m) ;
211  som_r_cheb (ti, nr, 1, 1, ksi, auxi) ;
212  coloc[i] = auxi[0]/
213  pow (-new_alpha*cos(M_PI*i/(nbre-1))+new_beta, power+l_quant);
214  }
215 
216  cfrcheb (dege, dege, coloc, dege, coloc) ;
217 
218  Tbl expansion (nbre) ;
219  expansion.set_etat_qcq() ;
220  for (int i=0 ; i<nbre ; i++) {
221  somme = 0 ;
222  for (int m=0 ; m<nbre ; m++)
223  somme += coloc[m]*tu(m, i) ;
224  expansion.set(i) = somme ;
225  }
226 
227  for (int i=0 ; i<nr ; i++) {
228  somme = 0 ;
229  for (int m=0 ; m<nbre ; m++)
230  somme += coef_u(m+l_quant, i)*expansion(m)*
231  pow(alpha_zec, m+l_quant)/
232  pow(new_alpha, m) ;
233  coef_zec.set(k, j, i) = somme ;
234  }
235  }
236  }
237 
238  va.set_etat_cf_qcq() ;
239  va.c_cf->set_etat_qcq() ;
240  va.c_cf->t[zone+1]->set_etat_qcq() ;
241 
242  for (int k=0 ; k<np+2 ; k++)
243  for (int j=0 ; j<nt ; j++)
244  for (int i=0 ; i<nr ; i++)
245  va.c_cf->set(zone+1, k, j, i) = coef_zec(k, j, i) ;
246 
247  set_dzpuis(power) ;
248  va.ylm_i() ;
249 
250  delete[] auxi ;
251  delete [] dege ;
252  delete [] ti ;
253  delete [] coloc ;
254 }
255 }
Bases of the spectral expansions.
Definition: base_val.h:322
const Map * mp
Reference mapping.
Definition: cmp.h:451
Valeur va
The numerical value of the Cmp
Definition: cmp.h:464
void set_dzpuis(int)
Set a value to dzpuis.
Definition: cmp.C:654
void raccord_externe(int puis, int nbre, int lmax)
Matching of the external domain with the outermost shell.
Affine radial mapping.
Definition: map.h:2027
const double * get_beta() const
Returns the pointer on the array beta.
Definition: map_af.C:481
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:477
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Matrix handling.
Definition: matrice.h:152
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:175
double & set(int j, int i)
Read/write of a particuliar element.
Definition: matrice.h:277
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
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
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:294
Basic array class.
Definition: tbl.h:161
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:372
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
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
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
Lorene prototypes.
Definition: app_hor.h:64