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 }
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Lorene::Base_val
Bases of the spectral expansions.
Definition: base_val.h:322
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Lorene::Valeur::set_etat_cf_qcq
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
Lorene::Matrice::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:175
Lorene::Cmp::raccord_c1_zec
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...
Definition: cmp_raccord_zec.C:84
Lorene::Mtbl_cf::t
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
Definition: mtbl_cf.h:205
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:372
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp
Definition: cmp.h:464
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Cmp::etat
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition: cmp.h:454
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Lorene::Mtbl_cf::set
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:294
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
Lorene::Matrice::inverse
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition: matrice.C:424
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Lorene::Matrice::set_band
void set_band(int up, int low) const
Calculate the band storage of *std.
Definition: matrice.C:364
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Lorene::Map_af
Affine radial mapping.
Definition: map.h:2027
Lorene::Cmp::dec2_dzpuis
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
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
Lorene::Matrice
Matrix handling.
Definition: matrice.h:152
Lorene::Base_val::get_base_r
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
Lorene::Cmp::mp
const Map * mp
Reference mapping.
Definition: cmp.h:451
Lorene::Valeur::base
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
Lorene::Matrice::set
double & set(int j, int i)
Read/write of a particuliar element.
Definition: matrice.h:277
Lorene::Cmp::dsdr
const Cmp & dsdr() const
Returns of *this .
Definition: cmp_deriv.C:84
Lorene::Map_af::get_alpha
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:477
Lorene::Cmp::set_dzpuis
void set_dzpuis(int)
Set a value to dzpuis.
Definition: cmp.C:654
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131
Lorene::Matrice::set_lu
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition: matrice.C:392