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 
23 
24 
25 /*
26  * $Id: cmp_raccord.C,v 1.5 2016/12/05 16:17:49 j_novak Exp $
27  * $Log: cmp_raccord.C,v $
28  * Revision 1.5 2016/12/05 16:17:49 j_novak
29  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
30  *
31  * Revision 1.4 2014/10/13 08:52:48 j_novak
32  * Lorene classes and functions now belong to the namespace Lorene.
33  *
34  * Revision 1.3 2014/10/06 15:13:03 j_novak
35  * Modified #include directives to use c++ syntax.
36  *
37  * Revision 1.2 2003/10/03 15:58:45 j_novak
38  * Cleaning of some headers
39  *
40  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
41  * LORENE
42  *
43  * Revision 2.1 2000/09/07 13:19:58 phil
44  * *** empty log message ***
45  *
46  * Revision 2.0 2000/06/06 12:18:27 phil
47  * *** empty log message ***
48  *
49  *
50  * $Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_raccord.C,v 1.5 2016/12/05 16:17:49 j_novak Exp $
51  *
52  */
53 
54 //standard
55 #include <cstdlib>
56 #include <cmath>
57 
58 // LORENE
59 #include "matrice.h"
60 #include "cmp.h"
61 #include "proto.h"
62 
63 
64 namespace Lorene {
65 Matrice matrice_raccord_pair (int cont, double alpha_kernel) {
66 
67  Matrice systeme (cont, cont) ;
68  systeme.set_etat_qcq() ;
69  for (int i=0 ; i<cont ; i++)
70  for (int j=0 ; j<cont ; j++)
71  systeme.set(i, j) = 0 ;
72 
73  double somme ;
74  for (int i=0 ; i<cont ; i++)
75  for (int k=0 ; k<cont ; k++)
76  if (k<= 2*i) {
77  somme = 1 ;
78  for (int boucle=0 ; boucle<k ; boucle++)
79  somme *= (4*i*i-boucle*boucle)/(2.*boucle+1.)/alpha_kernel ;
80  systeme.set(k, i) = somme ;
81  }
82  int inf = (cont%2 == 1) ? (cont-1)/2 : (cont-2)/2 ;
83  systeme.set_band (cont-1, inf) ;
84  systeme.set_lu() ;
85  return systeme ;
86 }
87 
88 Matrice matrice_raccord_impair (int cont, double alpha_kernel) {
89 
90  Matrice systeme (cont, cont) ;
91  systeme.set_etat_qcq() ;
92  for (int i=0 ; i<cont ; i++)
93  for (int j=0 ; j<cont ; j++)
94  systeme.set(i, j) = 0 ;
95 
96  double somme ;
97  for (int i=0 ; i<cont ; i++)
98  for (int k=0 ; k<cont ; k++)
99  if (k<= 2*i+1) {
100  somme = 1 ;
101  for (int boucle=0 ; boucle<k ; boucle++)
102  somme *= (pow(2*i+1, 2.)-boucle*boucle)/(2.*boucle+1.)/alpha_kernel ;
103  systeme.set(k, i) = somme ;
104  }
105  int inf = (cont%2 == 0) ? cont/2 : (cont-1)/2 ;
106  systeme.set_band (cont-1, inf) ;
107  systeme.set_lu() ;
108  return systeme ;
109 }
110 
111 
112 Tbl sec_membre_raccord (Tbl coef, int cont, double alpha_shell) {
113 
114  assert (coef.get_etat() != ETATNONDEF) ;
115  int nr = coef.get_dim(0) ;
116 
117  Tbl sec_membre(cont) ;
118  sec_membre.set_etat_qcq() ;
119  for (int i=0 ; i<cont ; i++)
120  sec_membre.set(i) = 0 ;
121 
122  double somme ;
123  for (int i=0 ; i<nr ; i++)
124  for (int k=0 ; k<cont ; k++)
125  if (k<= i) {
126  somme = 1 ;
127  for (int boucle=0 ; boucle<k ; boucle++)
128  somme *= (i*i-boucle*boucle)/(2.*boucle+1.)/alpha_shell ;
129  if ((i+k)%2 == 0)
130  sec_membre.set(k) += coef(i)*somme ;
131  else
132  sec_membre.set(k) -= coef(i)*somme ;
133  }
134 
135  return sec_membre ;
136 }
137 
138 
139 Tbl regularise (Tbl coef, int nr, int base_r) {
140 
141  assert ((base_r == R_CHEBI) || (base_r == R_CHEBP)) ;
142  assert (coef.get_etat() != ETATNONDEF) ;
143  int cont = coef.get_dim(0) ;
144  assert (nr >= cont) ;
145 
146  Tbl resu (nr) ;
147  resu.set_etat_qcq() ;
148 
149  double* x4coef = new double[nr] ;
150  for (int i=0 ; i<cont ; i++)
151  x4coef[i] = coef(i) ;
152  for (int i=cont ; i<nr ; i++)
153  x4coef[i] = 0 ;
154  double* x6coef = new double[nr] ;
155 
156  multx2_1d (nr, &x4coef, base_r) ;
157  multx2_1d (nr, &x4coef, base_r) ;
158  for (int i=0 ; i<nr ; i++)
159  x6coef[i] = x4coef[i] ;
160  multx2_1d (nr, &x6coef, base_r) ;
161 
162  for (int i=0 ; i<nr ; i++)
163  resu.set(i) = 3*x4coef[i]-2*x6coef[i] ;
164 
165  delete [] x4coef ;
166  delete [] x6coef ;
167 
168  return resu ;
169 }
170 
171 
172 
173 void Cmp::raccord (int aux) {
174  assert (etat != ETATNONDEF) ;
175 
176  assert (aux >=0) ;
177  int cont = aux+1 ;
178 
179  const Map_af* mapping = dynamic_cast<const Map_af*>(get_mp() ) ;
180 
181  if (mapping == 0x0) {
182  cout <<
183  "raccord : The mapping does not belong to the class Map_af !"
184  << endl ;
185  abort() ;
186  }
187 
188  assert (mapping->get_mg()->get_type_r(1) == FIN) ;
189  assert (mapping->get_mg()->get_type_r(0) == RARE) ;
190 
191  // On passe en Ylm et vire tout dans la zone interne...
192  va.coef() ;
193  va.ylm() ;
194  va.set_etat_cf_qcq() ;
195  va.c_cf->t[0]->annule_hard() ;
196 
197  // Confort :
198  int nz = mapping->get_mg()->get_nzone() ;
199  int nbrer_kernel = mapping->get_mg()->get_nr(0) ;
200  int nbrer_shell = mapping->get_mg()->get_nr(1) ;
201 
202  int nbret_kernel = mapping->get_mg()->get_nt(0) ;
203  int nbret_shell = mapping->get_mg()->get_nt(1) ;
204 
205  int nbrep_kernel = mapping->get_mg()->get_np(0) ;
206  int nbrep_shell = mapping->get_mg()->get_np(1) ;
207 
208  double alpha_kernel = mapping->get_alpha()[0] ;
209  double alpha_shell = mapping->get_alpha()[1] ;
210 
211  int base_r, m_quant, l_quant ;
212 
213  for (int k=0 ; k<nbrep_kernel+1 ; k++)
214  for (int j=0 ; j<nbret_kernel ; j++)
215  if (nullite_plm(j, nbret_kernel, k,nbrep_kernel, va.base) == 1)
216  if (nullite_plm(j, nbret_shell, k, nbrep_shell, va.base) == 1)
217  {
218  // calcul des nombres quantiques :
219  donne_lm(nz, 0, j, k, va.base, m_quant, l_quant, base_r) ;
220  assert ((base_r == R_CHEBP) || (base_r == R_CHEBI)) ;
221 
222  Matrice systeme(cont, cont) ;
223 
224  Tbl facteur (nbrer_kernel) ;
225  facteur.annule_hard() ;
226  for (int i=0 ; i<nbrer_shell ; i++)
227  if (i<nbrer_kernel)
228  facteur.set(i) = (*va.c_cf)(1, k, j, i) ;
229 
230  Tbl sec_membre (sec_membre_raccord (facteur, cont, alpha_shell)) ;
231 
232  if (base_r == R_CHEBP)
233  systeme = matrice_raccord_pair (cont, alpha_kernel) ;
234  else
235  systeme = matrice_raccord_impair (cont, alpha_kernel) ;
236 
237  Tbl soluce (systeme.inverse(sec_membre)) ;
238  Tbl regulier (nbrer_kernel) ;
239 
240  if (l_quant == 0)
241  for (int i=0 ; i<cont ; i++)
242  va.c_cf->set(0, k, j, i) = soluce(i) ;
243  else {
244  if (l_quant %2 == 0)
245  regulier = regularise (soluce, nbrer_kernel, R_CHEBP) ;
246  else
247  regulier = regularise (soluce, nbrer_kernel, R_CHEBI) ;
248 
249  for (int i=0 ; i<nbrer_kernel ; i++)
250  va.c_cf->set(0, k, j, i) = regulier(i) ;
251  }
252  }
253  va.ylm_i() ;
254 }
255 }
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:607
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:134
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
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:715
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:151
Lorene prototypes.
Definition: app_hor.h:67
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:304
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition: matrice.C:427
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:141
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition: cmp.h:454
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Matrix handling.
Definition: matrice.h:152
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:315
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:423
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Affine radial mapping.
Definition: map.h:2048
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Definition: cmp_raccord.C:173
Basic array class.
Definition: tbl.h:164
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:491
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:375
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s which contain the spectral coefficients in each domain...
Definition: mtbl_cf.h:215