LORENE
map_log_elliptic.C
1 
23 char map_log_elliptic_C[] = "$Header $" ;
24 
25 /*
26  * $Id: map_log_elliptic.C,v 1.6 2014/10/13 08:53:05 j_novak Exp $
27  * $Log: map_log_elliptic.C,v $
28  * Revision 1.6 2014/10/13 08:53:05 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.5 2014/10/06 15:13:13 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.4 2007/01/16 15:08:07 n_vasset
35  * New constructor, usn Scalar on mono-domain angular grid for boundary,
36  * for function sol_elliptic_boundary()
37  *
38  * Revision 1.3 2005/06/09 07:57:32 f_limousin
39  * Implement a new function sol_elliptic_boundary() and
40  * Vector::poisson_boundary(...) which solve the vectorial poisson
41  * equation (method 6) with an inner boundary condition.
42  *
43  * Revision 1.2 2004/08/24 09:14:42 p_grandclement
44  * Addition of some new operators, like Poisson in 2d... It now requieres the
45  * GSL library to work.
46  *
47  * Also, the way a variable change is stored by a Param_elliptic is changed and
48  * no longer uses Change_var but rather 2 Scalars. The codes using that feature
49  * will requiere some modification. (It should concern only the ones about monopoles)
50  *
51  * Revision 1.1 2004/06/22 08:49:58 p_grandclement
52  * Addition of everything needed for using the logarithmic mapping
53  *
54  * Revision 1.4 2004/03/17 15:58:48 p_grandclement
55  *
56  * $Header: /cvsroot/Lorene/C++/Source/Map/map_log_elliptic.C,v 1.6 2014/10/13 08:53:05 j_novak Exp $
57  *
58  */
59 
60 // Header C :
61 #include <cstdlib>
62 #include <cmath>
63 
64 // Headers Lorene :
65 #include "tbl.h"
66 #include "mtbl_cf.h"
67 #include "map.h"
68 #include "param_elliptic.h"
69 
70  //----------------------------------------------
71  // General elliptic solver
72  //----------------------------------------------
73 
74 namespace Lorene {
75 void Map_log::sol_elliptic(Param_elliptic& ope_var, const Scalar& source,
76  Scalar& pot) const {
77 
78  assert(source.get_etat() != ETATNONDEF) ;
79  assert(source.get_mp().get_mg() == mg) ;
80  assert(pot.get_mp().get_mg() == mg) ;
81 
82  assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
83  source.check_dzpuis(4)) ;
84  // Spherical harmonic expansion of the source
85  // ------------------------------------------
86 
87  const Valeur& sourva = source.get_spectral_va() ;
88 
89  if (sourva.get_etat() == ETATZERO) {
90  pot.set_etat_zero() ;
91  return ;
92  }
93 
94  // Spectral coefficients of the source
95  assert(sourva.get_etat() == ETATQCQ) ;
96 
97  Valeur rho(sourva.get_mg()) ;
98  sourva.coef() ;
99  rho = *(sourva.c_cf) ; // copy of the coefficients of the source
100 
101  rho.ylm() ; // spherical harmonic transforms
102 
103  // On met les bonnes bases dans le changement de variable
104  // de ope_var :
105  ope_var.var_F.set_spectral_va().coef() ;
106  ope_var.var_F.set_spectral_va().ylm() ;
107  ope_var.var_G.set_spectral_va().coef() ;
108  ope_var.var_G.set_spectral_va().ylm() ;
109 
110  // Call to the Mtbl_cf version
111  // ---------------------------
112  Mtbl_cf resu = elliptic_solver (ope_var, *(rho.c_cf)) ;
113  // Final result returned as a Scalar
114  // ---------------------------------
115 
116  pot.set_etat_zero() ; // to call Scalar::del_t().
117 
118  pot.set_etat_qcq() ;
119 
120  pot.set_spectral_va() = resu ;
121  pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
122 
123  pot.set_dzpuis(0) ;
124 }
125 
126 
127  //--------------------------------------------------
128  // General elliptic solver with boundary as Mtbl_cf
129  //--------------------------------------------------
130 
131 
133  Scalar& pot, const Mtbl_cf& bound,
134  double fact_dir, double fact_neu) const {
135 
136  assert(source.get_etat() != ETATNONDEF) ;
137  assert(source.get_mp().get_mg() == mg) ;
138  assert(pot.get_mp().get_mg() == mg) ;
139 
140  assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
141  source.check_dzpuis(4)) ;
142  // Spherical harmonic expansion of the source
143  // ------------------------------------------
144 
145  const Valeur& sourva = source.get_spectral_va() ;
146 
147  if (sourva.get_etat() == ETATZERO) {
148  pot.set_etat_zero() ;
149  return ;
150  }
151 
152  // Spectral coefficients of the source
153  assert(sourva.get_etat() == ETATQCQ) ;
154 
155  Valeur rho(sourva.get_mg()) ;
156  sourva.coef() ;
157  rho = *(sourva.c_cf) ; // copy of the coefficients of the source
158 
159  rho.ylm() ; // spherical harmonic transforms
160 
161  // On met les bonnes bases dans le changement de variable
162  // de ope_var :
163  ope_var.var_F.set_spectral_va().coef() ;
164  ope_var.var_F.set_spectral_va().ylm() ;
165  ope_var.var_G.set_spectral_va().coef() ;
166  ope_var.var_G.set_spectral_va().ylm() ;
167 
168  // Call to the Mtbl_cf version
169  // ---------------------------
170  Mtbl_cf resu = elliptic_solver_boundary (ope_var, *(rho.c_cf), bound,
171  fact_dir, fact_neu) ;
172  // Final result returned as a Scalar
173  // ---------------------------------
174 
175  pot.set_etat_zero() ; // to call Scalar::del_t().
176 
177  pot.set_etat_qcq() ;
178 
179  pot.set_spectral_va() = resu ;
180  pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
181 
182  pot.set_dzpuis(0) ;
183 }
184 
185 
186  //--------------------------------------------------
187  // General elliptic solver with boundary as Scalar
188  //--------------------------------------------------
189 
190 
192  Scalar& pot, const Scalar& bound,
193  double fact_dir, double fact_neu) const {
194 
195  assert(source.get_etat() != ETATNONDEF) ;
196  assert(source.get_mp().get_mg() == mg) ;
197  assert(pot.get_mp().get_mg() == mg) ;
198 
199  assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
200  source.check_dzpuis(4)) ;
201  // Spherical harmonic expansion of the source
202  // ------------------------------------------
203 
204  const Valeur& sourva = source.get_spectral_va() ;
205 
206  if (sourva.get_etat() == ETATZERO) {
207  pot.set_etat_zero() ;
208  return ;
209  }
210 
211  // Spectral coefficients of the source
212  assert(sourva.get_etat() == ETATQCQ) ;
213 
214  Valeur rho(sourva.get_mg()) ;
215  sourva.coef() ;
216  rho = *(sourva.c_cf) ; // copy of the coefficients of the source
217 
218  rho.ylm() ; // spherical harmonic transforms
219 
220  // On met les bonnes bases dans le changement de variable
221  // de ope_var :
222  ope_var.var_F.set_spectral_va().coef() ;
223  ope_var.var_F.set_spectral_va().ylm() ;
224  ope_var.var_G.set_spectral_va().coef() ;
225  ope_var.var_G.set_spectral_va().ylm() ;
226 
227  // Call to the Mtbl_cf version
228  // ---------------------------
229  Scalar bbound = bound;
230  bbound.set_spectral_va().ylm() ;
231  const Map& mapp = bbound.get_mp();
232 
233  const Mg3d& gri2d = *mapp.get_mg();
234 
235  assert(&gri2d == source.get_mp().get_mg()->get_angu_1dom()) ;
236 
237  Mtbl_cf bound2 (gri2d , bbound.get_spectral_base()) ;
238  bound2.annule_hard() ;
239 
240  if (bbound.get_etat() != ETATZERO){
241 
242  int nr = gri2d.get_nr(0) ;
243  int nt = gri2d.get_nt(0) ;
244  int np = gri2d.get_np(0) ;
245 
246  for(int k=0; k<np+2; k++)
247  for (int j=0; j<=nt-1; j++)
248  for(int xi=0; xi<= nr-1; xi++)
249  {
250  bound2.set(0, k , j , xi) = (*bbound.get_spectral_va().c_cf)(0, k, j, xi) ;
251  }
252  }
253 
254 
255  Mtbl_cf resu = elliptic_solver_boundary (ope_var, *(rho.c_cf), bound2,
256  fact_dir, fact_neu) ;
257  // Final result returned as a Scalar
258  // ---------------------------------
259 
260  pot.set_etat_zero() ; // to call Scalar::del_t().
261 
262  pot.set_etat_qcq() ;
263 
264  pot.set_spectral_va() = resu ;
265  pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
266 
267  pot.set_dzpuis(0) ;
268 }
269 
270 
271  //------------------------------------------------
272  // General elliptic solver with no zec
273  //-------------------------------------------------
274 
275 void Map_log::sol_elliptic_no_zec (Param_elliptic& ope_var, const Scalar& source,
276  Scalar& pot, double val) const {
277 
278  assert(source.get_etat() != ETATNONDEF) ;
279  assert(source.get_mp().get_mg() == mg) ;
280  assert(pot.get_mp().get_mg() == mg) ;
281 
282  assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
283  source.check_dzpuis(4)) ;
284  // Spherical harmonic expansion of the source
285  // ------------------------------------------
286 
287  const Valeur& sourva = source.get_spectral_va() ;
288 
289  if (sourva.get_etat() == ETATZERO) {
290  pot.set_etat_zero() ;
291  return ;
292  }
293 
294  // Spectral coefficients of the source
295  assert(sourva.get_etat() == ETATQCQ) ;
296 
297  Valeur rho(sourva.get_mg()) ;
298  sourva.coef() ;
299  rho = *(sourva.c_cf) ; // copy of the coefficients of the source
300 
301  rho.ylm() ; // spherical harmonic transforms
302 
303  // On met les bonnes bases dans le changement de variable
304  // de ope_var :
305  ope_var.var_F.set_spectral_va().coef() ;
306  ope_var.var_F.set_spectral_va().ylm() ;
307  ope_var.var_G.set_spectral_va().coef() ;
308  ope_var.var_G.set_spectral_va().ylm() ;
309 
310  // Call to the Mtbl_cf version
311  // ---------------------------
312  Mtbl_cf resu = elliptic_solver_no_zec (ope_var, *(rho.c_cf), val) ;
313  // Final result returned as a Scalar
314  // ---------------------------------
315 
316  pot.set_etat_zero() ; // to call Scalar::del_t().
317 
318  pot.set_etat_qcq() ;
319 
320  pot.set_spectral_va() = resu ;
321  pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
322 
323  pot.set_dzpuis(0) ;
324 }
325 
326 
327 }
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
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
Lorene::Valeur::get_etat
int get_etat() const
Returns the logical state.
Definition: valeur.h:726
Lorene::Mg3d
Multi-domain grid.
Definition: grilles.h:273
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::Param_elliptic::var_G
Scalar var_G
Multiplicative variable change that must be sphericaly symetric !
Definition: param_elliptic.h:145
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Map_log::sol_elliptic_boundary
void sol_elliptic_boundary(Param_elliptic &params, const Scalar &so, Scalar &uu, const Mtbl_cf &bound, double fact_dir, double fact_neu) const
General elliptic solver including inner boundary conditions.
Definition: map_log_elliptic.C:132
Lorene::Valeur::get_mg
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
Definition: valeur.h:729
Lorene::Scalar::check_dzpuis
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition: scalar.C:873
Lorene::Map_log::sol_elliptic_no_zec
void sol_elliptic_no_zec(Param_elliptic &params, const Scalar &so, Scalar &uu, double) const
General elliptic solver.
Definition: map_log_elliptic.C:275
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:861
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::Param_elliptic
This class contains the parameters needed to call the general elliptic solver.
Definition: param_elliptic.h:135
Lorene::Scalar::get_spectral_base
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va
Definition: scalar.h:1294
Lorene::Map_log::sol_elliptic
void sol_elliptic(Param_elliptic &params, const Scalar &so, Scalar &uu) const
General elliptic solver.
Definition: map_log_elliptic.C:75
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
Lorene::Scalar::get_spectral_va
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:601
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
Lorene::Param_elliptic::var_F
Scalar var_F
Additive variable change function.
Definition: param_elliptic.h:144
Lorene::Mtbl_cf::annule_hard
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition: mtbl_cf.C:312
Lorene::Mg3d::get_angu_1dom
const Mg3d * get_angu_1dom() const
Returns the pointer on the associated mono-domain angular grid.
Definition: mg3d.C:494
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::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:353
Lorene::Scalar::set_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:324
Lorene::Scalar::set_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:808
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
Lorene::Map
Base class for coordinate mappings.
Definition: map.h:670
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:186
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131
Lorene::Map::mg
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition: map.h:676