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
26char 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
59namespace Lorene {
60int 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
67
68 va.coef() ;
69 va.ylm() ;
70
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)*
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
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
223 va.ylm_i() ;
224
225 delete[] auxi ;
226 delete [] dege ;
227 delete [] ti ;
228 delete [] coloc ;
229}
230}
Bases of the spectral expansions.
Definition base_val.h:322
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
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
Matrix handling.
Definition matrice.h:152
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
friend Scalar cos(const Scalar &)
Cosine.
void raccord_externe(int puis, int nbre, int lmax)
Matching of the external domain with the outermost shell.
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
friend Scalar sqrt(const Scalar &)
Square root.
friend Scalar pow(const Scalar &, int)
Power .
Valeur va
The numerical value of the Scalar
Definition scalar.h:405
Basic array class.
Definition tbl.h:161
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
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Lorene prototypes.
Definition app_hor.h:64