LORENE
map_log_elliptic.C
1
23char 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
74namespace Lorene {
75void 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
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}
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.
void sol_elliptic(Param_elliptic &params, const Scalar &so, Scalar &uu) const
General elliptic solver.
void sol_elliptic_no_zec(Param_elliptic &params, const Scalar &so, Scalar &uu, double) const
General elliptic solver.
Base class for coordinate mappings.
Definition map.h:670
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition map.h:676
Multi-domain grid.
Definition grilles.h:273
const Mg3d * get_angu_1dom() const
Returns the pointer on the associated mono-domain angular grid.
Definition mg3d.C:494
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_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:294
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition mtbl_cf.C:312
This class contains the parameters needed to call the general elliptic solver.
Scalar var_G
Multiplicative variable change that must be sphericaly symetric !
Scalar var_F
Additive variable change function.
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
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
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition scalar.h:601
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va
Definition scalar.h:1294
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
Values and coefficients of a (real-value) function.
Definition valeur.h:287
int get_etat() const
Returns the logical state.
Definition valeur.h:726
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.
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
Definition valeur.h:729
void ylm_i()
Inverse of ylm()
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene prototypes.
Definition app_hor.h:64