LORENE
cmp_raccord.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
23char cmp_raccord_C[] = "$Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_raccord.C,v 1.4 2014/10/13 08:52:48 j_novak Exp $" ;
24
25/*
26 * $Id: cmp_raccord.C,v 1.4 2014/10/13 08:52:48 j_novak Exp $
27 * $Log: cmp_raccord.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:03 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.1 2000/09/07 13:19:58 phil
41 * *** empty log message ***
42 *
43 * Revision 2.0 2000/06/06 12:18:27 phil
44 * *** empty log message ***
45 *
46 *
47 * $Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_raccord.C,v 1.4 2014/10/13 08:52:48 j_novak Exp $
48 *
49 */
50
51//standard
52#include <cstdlib>
53#include <cmath>
54
55// LORENE
56#include "matrice.h"
57#include "cmp.h"
58#include "proto.h"
59
60
61namespace Lorene {
62Matrice matrice_raccord_pair (int cont, double alpha_kernel) {
63
64 Matrice systeme (cont, cont) ;
65 systeme.set_etat_qcq() ;
66 for (int i=0 ; i<cont ; i++)
67 for (int j=0 ; j<cont ; j++)
68 systeme.set(i, j) = 0 ;
69
70 double somme ;
71 for (int i=0 ; i<cont ; i++)
72 for (int k=0 ; k<cont ; k++)
73 if (k<= 2*i) {
74 somme = 1 ;
75 for (int boucle=0 ; boucle<k ; boucle++)
76 somme *= (4*i*i-boucle*boucle)/(2.*boucle+1.)/alpha_kernel ;
77 systeme.set(k, i) = somme ;
78 }
79 int inf = (cont%2 == 1) ? (cont-1)/2 : (cont-2)/2 ;
80 systeme.set_band (cont-1, inf) ;
81 systeme.set_lu() ;
82 return systeme ;
83}
84
85Matrice matrice_raccord_impair (int cont, double alpha_kernel) {
86
87 Matrice systeme (cont, cont) ;
88 systeme.set_etat_qcq() ;
89 for (int i=0 ; i<cont ; i++)
90 for (int j=0 ; j<cont ; j++)
91 systeme.set(i, j) = 0 ;
92
93 double somme ;
94 for (int i=0 ; i<cont ; i++)
95 for (int k=0 ; k<cont ; k++)
96 if (k<= 2*i+1) {
97 somme = 1 ;
98 for (int boucle=0 ; boucle<k ; boucle++)
99 somme *= (pow(2*i+1, 2.)-boucle*boucle)/(2.*boucle+1.)/alpha_kernel ;
100 systeme.set(k, i) = somme ;
101 }
102 int inf = (cont%2 == 0) ? cont/2 : (cont-1)/2 ;
103 systeme.set_band (cont-1, inf) ;
104 systeme.set_lu() ;
105 return systeme ;
106}
107
108
109Tbl sec_membre_raccord (Tbl coef, int cont, double alpha_shell) {
110
111 assert (coef.get_etat() != ETATNONDEF) ;
112 int nr = coef.get_dim(0) ;
113
114 Tbl sec_membre(cont) ;
115 sec_membre.set_etat_qcq() ;
116 for (int i=0 ; i<cont ; i++)
117 sec_membre.set(i) = 0 ;
118
119 double somme ;
120 for (int i=0 ; i<nr ; i++)
121 for (int k=0 ; k<cont ; k++)
122 if (k<= i) {
123 somme = 1 ;
124 for (int boucle=0 ; boucle<k ; boucle++)
125 somme *= (i*i-boucle*boucle)/(2.*boucle+1.)/alpha_shell ;
126 if ((i+k)%2 == 0)
127 sec_membre.set(k) += coef(i)*somme ;
128 else
129 sec_membre.set(k) -= coef(i)*somme ;
130 }
131
132 return sec_membre ;
133}
134
135
136Tbl regularise (Tbl coef, int nr, int base_r) {
137
138 assert ((base_r == R_CHEBI) || (base_r == R_CHEBP)) ;
139 assert (coef.get_etat() != ETATNONDEF) ;
140 int cont = coef.get_dim(0) ;
141 assert (nr >= cont) ;
142
143 Tbl resu (nr) ;
144 resu.set_etat_qcq() ;
145
146 double* x4coef = new double[nr] ;
147 for (int i=0 ; i<cont ; i++)
148 x4coef[i] = coef(i) ;
149 for (int i=cont ; i<nr ; i++)
150 x4coef[i] = 0 ;
151 double* x6coef = new double[nr] ;
152
153 multx2_1d (nr, &x4coef, base_r) ;
154 multx2_1d (nr, &x4coef, base_r) ;
155 for (int i=0 ; i<nr ; i++)
156 x6coef[i] = x4coef[i] ;
157 multx2_1d (nr, &x6coef, base_r) ;
158
159 for (int i=0 ; i<nr ; i++)
160 resu.set(i) = 3*x4coef[i]-2*x6coef[i] ;
161
162 delete [] x4coef ;
163 delete [] x6coef ;
164
165 return resu ;
166}
167
168
169
170void Cmp::raccord (int aux) {
171 assert (etat != ETATNONDEF) ;
172
173 assert (aux >=0) ;
174 int cont = aux+1 ;
175
176 const Map_af* mapping = dynamic_cast<const Map_af*>(get_mp() ) ;
177
178 if (mapping == 0x0) {
179 cout <<
180 "raccord : The mapping does not belong to the class Map_af !"
181 << endl ;
182 abort() ;
183 }
184
185 assert (mapping->get_mg()->get_type_r(1) == FIN) ;
186 assert (mapping->get_mg()->get_type_r(0) == RARE) ;
187
188 // On passe en Ylm et vire tout dans la zone interne...
189 va.coef() ;
190 va.ylm() ;
192 va.c_cf->t[0]->annule_hard() ;
193
194 // Confort :
195 int nz = mapping->get_mg()->get_nzone() ;
196 int nbrer_kernel = mapping->get_mg()->get_nr(0) ;
197 int nbrer_shell = mapping->get_mg()->get_nr(1) ;
198
199 int nbret_kernel = mapping->get_mg()->get_nt(0) ;
200 int nbret_shell = mapping->get_mg()->get_nt(1) ;
201
202 int nbrep_kernel = mapping->get_mg()->get_np(0) ;
203 int nbrep_shell = mapping->get_mg()->get_np(1) ;
204
205 double alpha_kernel = mapping->get_alpha()[0] ;
206 double alpha_shell = mapping->get_alpha()[1] ;
207
208 int base_r, m_quant, l_quant ;
209
210 for (int k=0 ; k<nbrep_kernel+1 ; k++)
211 for (int j=0 ; j<nbret_kernel ; j++)
212 if (nullite_plm(j, nbret_kernel, k,nbrep_kernel, va.base) == 1)
213 if (nullite_plm(j, nbret_shell, k, nbrep_shell, va.base) == 1)
214 {
215 // calcul des nombres quantiques :
216 donne_lm(nz, 0, j, k, va.base, m_quant, l_quant, base_r) ;
217 assert ((base_r == R_CHEBP) || (base_r == R_CHEBI)) ;
218
219 Matrice systeme(cont, cont) ;
220
221 Tbl facteur (nbrer_kernel) ;
222 facteur.annule_hard() ;
223 for (int i=0 ; i<nbrer_shell ; i++)
224 if (i<nbrer_kernel)
225 facteur.set(i) = (*va.c_cf)(1, k, j, i) ;
226
227 Tbl sec_membre (sec_membre_raccord (facteur, cont, alpha_shell)) ;
228
229 if (base_r == R_CHEBP)
230 systeme = matrice_raccord_pair (cont, alpha_kernel) ;
231 else
232 systeme = matrice_raccord_impair (cont, alpha_kernel) ;
233
234 Tbl soluce (systeme.inverse(sec_membre)) ;
235 Tbl regulier (nbrer_kernel) ;
236
237 if (l_quant == 0)
238 for (int i=0 ; i<cont ; i++)
239 va.c_cf->set(0, k, j, i) = soluce(i) ;
240 else {
241 if (l_quant %2 == 0)
242 regulier = regularise (soluce, nbrer_kernel, R_CHEBP) ;
243 else
244 regulier = regularise (soluce, nbrer_kernel, R_CHEBI) ;
245
246 for (int i=0 ; i<nbrer_kernel ; i++)
247 va.c_cf->set(0, k, j, i) = regulier(i) ;
248 }
249 }
250 va.ylm_i() ;
251}
252}
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition cmp.h:454
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Affine radial mapping.
Definition map.h:2027
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
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition matrice.C:424
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
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:474
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
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_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 pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define R_CHEBP
base de Cheb. paire (rare) seulement
Lorene prototypes.
Definition app_hor.h:64