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