LORENE
cmp_raccord_zec.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 char cmp_raccord_zec_C[] = "$Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_raccord_zec.C,v 1.4 2014/10/13 08:52:48 j_novak Exp $" ;
24 
25 /*
26  * $Id: cmp_raccord_zec.C,v 1.4 2014/10/13 08:52:48 j_novak Exp $
27  * $Log: cmp_raccord_zec.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.7 2001/03/30 13:38:32 phil
41  * *** empty log message ***
42  *
43  * Revision 2.6 2001/03/22 10:25:01 phil
44  * changement complet : cas plus general
45  *
46  * Revision 2.5 2001/02/08 14:21:32 phil
47  * correction de raccord_zec.C (on prend en compte le dernier coef ...)
48  *
49  * Revision 2.4 2001/01/02 11:25:37 phil
50  * *** empty log message ***
51  *
52  * Revision 2.3 2000/12/13 14:59:18 phil
53  * *** empty log message ***
54  *
55  * Revision 2.2 2000/12/13 14:49:54 phil
56  * changement nom variable appel
57  * /
58  *
59  * Revision 2.1 2000/12/13 14:12:29 phil
60  * correction bugs
61  *
62  * Revision 2.0 2000/12/13 14:09:31 phil
63  * *** empty log message ***
64  *
65  *
66  * $Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_raccord_zec.C,v 1.4 2014/10/13 08:52:48 j_novak Exp $
67  *
68  */
69 
70 //standard
71 #include <cstdlib>
72 #include <cmath>
73 
74 // LORENE
75 #include "matrice.h"
76 #include "cmp.h"
77 #include "proto.h"
78 
79 // Fait le raccord C1 dans la zec ...
80 namespace Lorene {
81 // Suppose (pour le moment, le meme nbre de points sur les angles ...)
82 // et que la zone precedente est une coquille
83 
84 void Cmp::raccord_c1_zec (int puis, int nbre, int lmax) {
85 
86  assert (nbre>0) ;
87  assert (etat != ETATNONDEF) ;
88  if (etat == ETATZERO)
89  return ;
90 
91  // Le mapping doit etre affine :
92  const Map_af* map = dynamic_cast<const Map_af*>(mp) ;
93  if (map == 0x0) {
94  cout << "Le mapping doit etre affine" << endl ;
95  abort() ;
96  }
97 
98  int nz = map->get_mg()->get_nzone() ;
99  int nr = map->get_mg()->get_nr (nz-1) ;
100  int nt = map->get_mg()->get_nt (nz-1) ;
101  int np = map->get_mg()->get_np (nz-1) ;
102 
103  double alpha = map->get_alpha()[nz-1] ;
104  double r_cont = -1./2./alpha ; //Rayon de debut de la zec.
105 
106  // On calcul les coefficients des puissances de 1./r
107  Tbl coef (nbre+2*lmax, nr) ;
108  coef.set_etat_qcq() ;
109 
110  int* deg = new int[3] ;
111  deg[0] = 1 ; deg[1] = 1 ; deg[2] = nr ;
112  double* auxi = new double[nr] ;
113 
114  for (int conte=0 ; conte<nbre+2*lmax ; conte++) {
115  for (int i=0 ; i<nr ; i++)
116  auxi[i] = pow(-1-cos(M_PI*i/(nr-1)), (conte+puis)) ;
117 
118  cfrcheb(deg, deg, auxi, deg, auxi) ;
119  for (int i=0 ; i<nr ; i++)
120  coef.set(conte, i) = auxi[i]*pow (alpha, conte+puis) ;
121  }
122 
123  delete[] deg ;
124  // Maintenant on va calculer les valeurs de la ieme derivee :
125  Tbl valeurs (nbre, nt, np+1) ;
126  valeurs.set_etat_qcq() ;
127 
128  Cmp courant (*this) ;
129  double* res_val = new double[1] ;
130 
131  for (int conte=0 ; conte<nbre ; conte++) {
132 
133  courant.va.coef() ;
134  courant.va.ylm() ;
135  courant.va.c_cf->t[nz-1]->annule_hard() ;
136 
137  int base_r = courant.va.base.get_base_r(nz-2) ;
138  for (int k=0 ; k<np+1 ; k++)
139  for (int j=0 ; j<nt ; j++)
140  if (nullite_plm(j, nt, k, np, courant.va.base) == 1) {
141 
142  for (int i=0 ; i<nr ; i++)
143  auxi[i] = (*courant.va.c_cf)(nz-2, k, j, i) ;
144 
145  switch (base_r) {
146  case R_CHEB :
147  som_r_cheb (auxi, nr, 1, 1, 1, res_val) ;
148  break ;
149  default :
150  cout << "Cas non prevu dans raccord_zec" << endl ;
151  abort() ;
152  break ;
153  }
154  valeurs.set(conte, k, j) = res_val[0] ;
155  }
156  Cmp copie (courant) ;
157  copie.dec2_dzpuis() ;
158  courant = copie.dsdr() ;
159  }
160 
161  delete [] auxi ;
162  delete [] res_val ;
163 
164  // On boucle sur les harmoniques : construction de la matrice
165  // et du second membre
166  va.coef() ;
167  va.ylm() ;
168  va.c_cf->t[nz-1]->annule_hard() ;
169  va.set_etat_cf_qcq() ;
170 
171  const Base_val& base = va.base ;
172  int base_r, l_quant, m_quant ;
173  for (int k=0 ; k<np+1 ; k++)
174  for (int j=0 ; j<nt ; j++)
175  if (nullite_plm(j, nt, k, np, va.base) == 1) {
176 
177  donne_lm (nz, nz-1, j, k, base, m_quant, l_quant, base_r) ;
178 
179  if (l_quant<=lmax) {
180 
181  Matrice systeme (nbre, nbre) ;
182  systeme.set_etat_qcq() ;
183 
184  for (int col=0 ; col<nbre ; col++)
185  for (int lig=0 ; lig<nbre ; lig++) {
186 
187  int facteur = (lig%2==0) ? 1 : -1 ;
188  for (int conte=0 ; conte<lig ; conte++)
189  facteur *= puis+col+conte+2*l_quant ;
190  systeme.set(lig, col) = facteur/pow(r_cont, puis+col+lig+2*l_quant) ;
191  }
192 
193  systeme.set_band(nbre, nbre) ;
194  systeme.set_lu() ;
195 
196  Tbl sec_membre (nbre) ;
197  sec_membre.set_etat_qcq() ;
198  for (int conte=0 ; conte<nbre ; conte++)
199  sec_membre.set(conte) = valeurs(conte, k, j) ;
200 
201  Tbl inv (systeme.inverse(sec_membre)) ;
202 
203  for (int conte=0 ; conte<nbre ; conte++)
204  for (int i=0 ; i<nr ; i++)
205  va.c_cf->set(nz-1, k, j, i)+=
206  inv(conte)*coef(conte+2*l_quant, i) ;
207  }
208  else for (int i=0 ; i<nr ; i++)
209  va.c_cf->set(nz-1, k, j, i)
210  = 0 ;
211  }
212 
213  va.ylm_i() ;
214  set_dzpuis (0) ;
215 }
216 }
Bases of the spectral expansions.
Definition: base_val.h:322
int get_base_r(int l) const
Returns the expansion basis for r ( ) functions in the domain of index l (e.g.
Definition: base_val.h:400
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
const Map * mp
Reference mapping.
Definition: cmp.h:451
void raccord_c1_zec(int puis, int nbre, int lmax)
Performs the matching of the external domain with respect to the last shell using function like wit...
Valeur va
The numerical value of the Cmp
Definition: cmp.h:464
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition: cmp.h:454
void set_dzpuis(int)
Set a value to dzpuis.
Definition: cmp.C:654
void dec2_dzpuis()
Decreases by 2 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
Definition: cmp_r_manip.C:180
const Cmp & dsdr() const
Returns of *this .
Definition: cmp_deriv.C:84
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
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
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition: matrice.C:424
void set_band(int up, int low) const
Calculate the band storage of *std.
Definition: matrice.C:364
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition: matrice.C:392
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
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 pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
Lorene prototypes.
Definition: app_hor.h:64