LORENE
cmp_manip.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_manip.C,v 1.7 2016/12/05 16:17:48 j_novak Exp $
27  * $Log: cmp_manip.C,v $
28  * Revision 1.7 2016/12/05 16:17:48 j_novak
29  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
30  *
31  * Revision 1.6 2014/10/13 08:52:47 j_novak
32  * Lorene classes and functions now belong to the namespace Lorene.
33  *
34  * Revision 1.5 2014/10/06 15:13:03 j_novak
35  * Modified #include directives to use c++ syntax.
36  *
37  * Revision 1.4 2008/08/19 06:41:59 j_novak
38  * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
39  * cast-type operations, and constant strings that must be defined as const char*
40  *
41  * Revision 1.3 2003/10/23 09:41:27 p_grandclement
42  * small modif of set_val_hor (one can work at the origin now)
43  *
44  * Revision 1.2 2003/10/03 15:58:44 j_novak
45  * Cleaning of some headers
46  *
47  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
48  * LORENE
49  *
50  * Revision 2.2 2001/05/25 09:29:58 phil
51  * ajout de filtre_phi
52  *
53  * Revision 2.1 2001/02/12 18:08:51 phil
54  * ajout de fixe_decroissance
55  *
56  * Revision 2.0 2000/10/19 09:23:37 phil
57  * *** empty log message ***
58  *
59  *
60  * $Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_manip.C,v 1.7 2016/12/05 16:17:48 j_novak Exp $
61  *
62  */
63 
64 //standard
65 #include <cstdlib>
66 #include <cmath>
67 
68 // Lorene
69 #include "cmp.h"
70 #include "proto.h"
71 
72 /*
73  * Annule les n derniers coefficients en r dans la derniere zone
74  */
75 
76 namespace Lorene {
77 void Cmp::filtre (int n) {
78 
79  assert (etat != ETATNONDEF) ;
80  if (etat == ETATZERO)
81  return ;
82 
83  int nz = mp->get_mg()->get_nzone() ;
84  int np = mp->get_mg()->get_np(nz-1) ;
85  int nt = mp->get_mg()->get_nt(nz-1) ;
86  int nr = mp->get_mg()->get_nr(nz-1) ;
87 
88  del_deriv() ;
89 
90  va.coef() ;
92 
93  for (int k=0 ; k<np+1 ; k++)
94  if (k!=1)
95  for (int j=0 ; j<nt ; j++)
96  for (int i=nr-1 ; i>nr-1-n ; i--)
97  va.c_cf->set(nz-1, k, j, i) = 0 ;
98 }
99 
100 /*
101  * Annule les n derniers coefficients en phi dans zone nz
102  */
103 
104 void Cmp::filtre_phi (int n, int nz) {
105  assert (etat != ETATNONDEF) ;
106  if (etat == ETATZERO)
107  return ;
108 
109  del_deriv() ;
110 
111  va.coef() ;
112  va.set_etat_cf_qcq() ;
113  int np = mp->get_mg()->get_np(nz) ;
114  int nt = mp->get_mg()->get_nt(nz) ;
115  int nr = mp->get_mg()->get_nr(nz) ;
116 
117  for (int k=np+1-n ; k<np+1 ; k++)
118  for (int j=0 ; j<nt ; j++)
119  for (int i=0 ; i<nr ; i++)
120  va.c_cf->set(nz, k, j, i) = 0 ;
121 }
122 
123 /*
124  * Fixe la valeur a l'infini (si la derniere zone est compactifiee)
125  * d'un Cmp a val
126  * Utile quand on a affaire a des nan0x10000000
127  */
128 
129 void Cmp::set_val_inf (double val) {
130 
131  assert (etat != ETATNONDEF) ;
132  if (etat == ETATZERO) {
133  if (val == 0)
134  return ;
135  else
136  annule_hard() ;
137  }
138  del_deriv() ;
139 
140  int nz = mp->get_mg()->get_nzone() ;
141 
142  // On verifie la compactification
143  assert (mp->get_mg()->get_type_r(nz-1) == UNSURR) ;
144 
145  int nr = mp->get_mg()->get_nr(nz-1) ;
146  int nt = mp->get_mg()->get_nt(nz-1) ;
147  int np = mp->get_mg()->get_np(nz-1) ;
148 
149  va.coef_i() ;
150  va.set_etat_c_qcq() ;
151 
152  for (int k=0 ; k<np ; k++)
153  for (int j=0 ; j<nt ; j++)
154  va.set(nz-1, k, j, nr-1) = val ;
155 }
156 
157 /*
158  * Fixe la valeur d'un Cmp a val, sur la frontiere interne de la coquille zone.
159  * Utile quand on a affaire a des nan0x10000000
160  */
161 
162 void Cmp::set_val_hor (double val, int zone) {
163 
164  assert (etat != ETATNONDEF) ;
165  if (etat == ETATZERO) {
166  if (val == 0)
167  return ;
168  else
169  annule_hard() ;
170  }
171  assert (zone < mp->get_mg()->get_nzone()) ;
172  del_deriv() ;
173 
174  int nt = mp->get_mg()->get_nt(zone) ;
175  int np = mp->get_mg()->get_np(zone) ;
176 
177  va.coef_i() ;
178  va.set_etat_c_qcq() ;
179 
180  for (int k=0 ; k<np ; k++)
181  for (int j=0 ; j<nt ; j++)
182  va.set(zone, k, j, 0) = val ;
183 }
184 
185 /*
186  * Permet de fixer la decroissance du cmp a l infini en viurant les
187  * termes en 1/r^n
188  */
189 void Cmp::fixe_decroissance (int puis) {
190 
191  if (puis<dzpuis)
192  return ;
193  else {
194 
195  int nbre = puis-dzpuis ;
196 
197  // le confort avant tout ! (c'est bien le confort ...)
198  int nz = mp->get_mg()->get_nzone() ;
199  int np = mp->get_mg()->get_np(nz-1) ;
200  int nt = mp->get_mg()->get_nt(nz-1) ;
201  int nr = mp->get_mg()->get_nr(nz-1) ;
202 
203  const Map_af* map = dynamic_cast<const Map_af*>(mp) ;
204  if (map == 0x0) {
205  cout << "Le mapping doit etre affine" << endl ;
206  abort() ;
207  }
208 
209  double alpha = map->get_alpha()[nz-1] ;
210 
211  Cmp courant (*this) ;
212 
213  va.coef() ;
214  va.set_etat_cf_qcq() ;
215 
216  for (int conte=0 ; conte<nbre ; conte++) {
217 
218  int base_r = courant.va.base.get_base_r(nz-1) ;
219 
220  courant.va.coef() ;
221 
222  // On calcul les coefficients de 1/r^conte
223  double* coloc = new double [nr] ;
224  int * deg = new int[3] ;
225  deg[0] = 1 ;
226  deg[1] = 1 ;
227  deg[2] = nr ;
228 
229  for (int i=0 ; i<nr ; i++)
230  coloc[i] =pow(alpha, double(conte))*
231  pow(-1-cos(M_PI*i/(nr-1)), double(conte)) ;
232 
233  cfrcheb(deg, deg, coloc, deg, coloc) ;
234 
235  for (int k=0 ; k<np+1 ; k++)
236  if (k != 1)
237  for (int j=0 ; j<nt ; j++) {
238 
239  // On doit determiner le coefficient du truc courant :
240  double* coef = new double [nr] ;
241  double* auxi = new double[1] ;
242  for (int i=0 ; i<nr ; i++)
243  coef[i] = (*courant.va.c_cf)(nz-1, k, j, i) ;
244  switch (base_r) {
245  case R_CHEBU :
246  som_r_chebu (coef, nr, 1, 1, 1, auxi) ;
247  break ;
248  default :
249  som_r_pas_prevu (coef, nr, 1, 1, 1, auxi) ;
250  break ;
251  }
252 
253  // On modifie le cmp courant :
254  courant.va.coef() ;
255  courant.va.set_etat_cf_qcq() ;
256  courant.va.c_cf->set(nz-1, k, j, 0) -= *auxi ;
257 
258  for (int i=0 ; i<nr ; i++)
259  this->va.c_cf->set(nz-1, k, j, i) -= *auxi * coloc[i] ;
260 
261 
262  delete [] coef ;
263  delete [] auxi ;
264  }
265  delete [] coloc ;
266  delete [] deg ;
267 
268  courant.mult_r_zec() ;
269  }
270  }
271 }
272 }
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:607
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
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
void coef_i() const
Computes the physical value of *this.
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition: cmp.h:454
void mult_r_zec()
Multiplication by r in the external compactified domain (ZEC)
Definition: cmp_r_manip.C:106
void del_deriv()
Logical destructor of the derivatives.
Definition: cmp.C:268
void annule_hard()
Sets the Cmp to zero in a hard way.
Definition: cmp.C:341
void filtre_phi(int n, int zone)
Sets the n lasts coefficients in to 0 in the domain zone .
Definition: cmp_manip.C:104
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:97
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:403
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Definition: cmp_manip.C:77
int dzpuis
Power of r by which the quantity represented by this must be divided in the external compactified zo...
Definition: cmp.h:461
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:315
void set_val_inf(double val)
Sets the value of the Cmp to val at infinity.
Definition: cmp_manip.C:129
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
void fixe_decroissance(int puis)
Substracts all the components behaving like in the external domain, with n strictly lower than puis ...
Definition: cmp_manip.C:189
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
void set_val_hor(double val, int zone)
Sets the value of the Cmp to val on the inner boudary of the shell number zone .This is usefull for d...
Definition: cmp_manip.C:162
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
const Map * mp
Reference mapping.
Definition: cmp.h:451
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c )...
Definition: valeur.C:704
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:491
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition: valeur.h:373