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
23char 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
67namespace Lorene {
68int 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
91void 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
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 & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:294
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
Basic array class.
Definition tbl.h:161
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 & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
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.
void ylm_i()
Inverse of ylm()
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