LORENE
vector_poisson.C
1 /*
2  * Methods for solving vector Poisson equation.
3  *
4  * (see file vector.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2003 Eric Gourgoulhon & Jerome Novak
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 char vector_poisson_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/vector_poisson.C,v 1.24 2014/10/13 08:53:45 j_novak Exp $" ;
29 
30 /*
31  * $Id: vector_poisson.C,v 1.24 2014/10/13 08:53:45 j_novak Exp $
32  * $Log: vector_poisson.C,v $
33  * Revision 1.24 2014/10/13 08:53:45 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.23 2014/10/06 15:13:21 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.22 2005/02/14 13:01:50 j_novak
40  * p_eta and p_mu are members of the class Vector. Most of associated functions
41  * have been moved from the class Vector_divfree to the class Vector.
42  *
43  * Revision 1.21 2005/01/10 15:36:09 j_novak
44  * In method 5: added transformation back from the Ylm base.
45  *
46  * Revision 1.20 2004/12/23 10:23:06 j_novak
47  * Improved method 5 in the case when some components are zero.
48  * Changed Vector_divfree::poisson to deduce eta from the equation.
49  *
50  * Revision 1.19 2004/08/24 09:14:50 p_grandclement
51  * Addition of some new operators, like Poisson in 2d... It now requieres the
52  * GSL library to work.
53  *
54  * Also, the way a variable change is stored by a Param_elliptic is changed and
55  * no longer uses Change_var but rather 2 Scalars. The codes using that feature
56  * will requiere some modification. (It should concern only the ones about monopoles)
57  *
58  * Revision 1.18 2004/07/27 09:40:05 j_novak
59  * Yet another method for solving vector Poisson eq. with spherical coordinates.
60  *
61  * Revision 1.17 2004/06/14 15:15:58 j_novak
62  * New method (No.4) to solve the vector Poisson eq. The output is continuous
63  * for all (spherical) components.
64  *
65  * Revision 1.16 2004/05/26 07:30:59 j_novak
66  * Version 1.15 was not good.
67  *
68  * Revision 1.15 2004/05/25 15:15:47 f_limousin
69  * Function Vector::poisson with parameters now returns a Vector (the
70  * result) instead of void.
71  *
72  * Revision 1.12 2004/03/26 17:05:24 j_novak
73  * Added new method n.3 using Tenseur::poisson_vect_oohara
74  *
75  * Revision 1.11 2004/03/11 08:48:45 f_limousin
76  * Implement method Vector::poisson with parameters, only with method
77  * 2 yet.
78  *
79  * Revision 1.10 2004/03/10 16:38:38 e_gourgoulhon
80  * Modified the prototype of poisson with param. to let it
81  * agree with declaration in vector.h.
82  *
83  * Revision 1.9 2004/03/03 09:07:03 j_novak
84  * In Vector::poisson(double, int), the flat metric is taken from the mapping.
85  *
86  * Revision 1.8 2004/02/24 17:00:25 j_novak
87  * Added a forgotten term.
88  *
89  * Revision 1.7 2004/02/24 09:46:20 j_novak
90  * Correction to cope with SGI compiler's warnings.
91  *
92  * Revision 1.6 2004/02/22 15:47:46 j_novak
93  * Added 2 more methods to solve the vector poisson equation. Method 1 is not
94  * tested yet.
95  *
96  * Revision 1.5 2004/02/20 10:53:41 j_novak
97  * Minor modifs.
98  *
99  * Revision 1.4 2004/02/16 17:40:14 j_novak
100  * Added a version of poisson with the flat metric as argument (avoids
101  * unnecessary calculations by decompose_div)
102  *
103  * Revision 1.3 2003/10/29 11:04:34 e_gourgoulhon
104  * dec2_dpzuis() replaced by dec_dzpuis(2).
105  * inc2_dpzuis() replaced by inc_dzpuis(2).
106  *
107  * Revision 1.2 2003/10/22 13:08:06 j_novak
108  * Better handling of dzpuis flags
109  *
110  * Revision 1.1 2003/10/20 15:15:42 j_novak
111  * New method Vector::poisson().
112  *
113  *
114  * $Headers: $
115  *
116  */
117 
118 //C headers
119 #include <cstdlib>
120 #include <cmath>
121 
122 // Lorene headers
123 #include "metric.h"
124 #include "tenseur.h"
125 #include "param.h"
126 #include "param_elliptic.h"
127 #include "proto.h"
128 
129 namespace Lorene {
130 Vector Vector::poisson(double lambda, const Metric_flat& met_f, int method)
131  const {
132 
133  bool nullite = true ;
134  for (int i=0; i<3; i++) {
135  assert(cmp[i]->check_dzpuis(4)) ;
136  if (cmp[i]->get_etat() != ETATZERO) nullite = false ;
137  }
138  assert ((method>=0) && (method<7)) ;
139 
140  Vector resu(*mp, CON, triad) ;
141  if (nullite)
142  resu.set_etat_zero() ;
143  else {
144 
145  switch (method) {
146 
147  case 0 : {
148 
149  Scalar poten(*mp) ;
150  if (fabs(lambda+1) < 1.e-8)
151  poten.set_etat_zero() ;
152  else {
153  poten = (potential(met_f) / (lambda + 1)).poisson() ;
154  }
155 
156  Vector grad = poten.derive_con(met_f) ;
157  grad.dec_dzpuis(2) ;
158 
159  return ( div_free(met_f).poisson() + grad) ;
160  break ;
161  }
162 
163  case 1 : {
164 
165  Scalar divf(*mp) ;
166  if (fabs(lambda+1) < 1.e-8)
167  divf.set_etat_zero() ;
168  else {
169  divf = (potential(met_f) / (lambda + 1)) ;
170  }
171 
172  int nz = mp->get_mg()->get_nzone() ;
173 
174  //-----------------------------------
175  // Removal of the l=0 part of div(F)
176  //-----------------------------------
177  Scalar div_lnot0 = divf ;
178  div_lnot0.div_r_dzpuis(4) ;
179  Scalar source_r(*mp) ;
180  Valeur& va_div = div_lnot0.set_spectral_va() ;
181  if (div_lnot0.get_etat() != ETATZERO) {
182  va_div.coef() ;
183  va_div.ylm() ;
184  for (int lz=0; lz<nz; lz++) {
185  int np = mp->get_mg()->get_np(lz) ;
186  int nt = mp->get_mg()->get_nt(lz) ;
187  int nr = mp->get_mg()->get_nr(lz) ;
188  if (va_div.c_cf->operator()(lz).get_etat() != ETATZERO)
189  for (int k=0; k<np+1; k++)
190  for (int j=0; j<nt; j++) {
191  int l_q, m_q, base_r ;
192  if (nullite_plm(j, nt, k, np, va_div.base) == 1) {
193  donne_lm(nz, lz, j, k, va_div.base, m_q, l_q, base_r) ;
194  if (l_q == 0)
195  for (int i=0; i<nr; i++)
196  va_div.c_cf->set(lz, k, j, i) = 0. ;
197  }
198  }
199  }
200  source_r.set_etat_qcq() ;
201  source_r.set_spectral_va() = 2*(*va_div.c_cf) ; //2*div(F)
202  source_r.set_spectral_va().ylm_i() ;
203  source_r.set_dzpuis(4) ;
204  }
205  else
206  source_r.set_etat_zero() ;
207 
208  //------------------------
209  // Other source terms ....
210  //------------------------
211  source_r += *(cmp[0]) - lambda*divf.dsdr() ;
212  Scalar f_r(*mp) ;
213  if (source_r.get_etat() != ETATZERO) {
214 
215  //------------------------
216  // The elliptic operator
217  //------------------------
218  Param_elliptic param_fr(source_r) ;
219  for (int lz=0; lz<nz; lz++)
220  param_fr.set_poisson_vect_r(lz) ;
221 
222  f_r = source_r.sol_elliptic(param_fr) ;
223  }
224  else
225  f_r.set_etat_zero() ;
226 
227  divf.dec_dzpuis(1) ;
228  Scalar source_eta = divf - f_r.dsdr() ;
229  source_eta.mult_r_dzpuis(0) ;
230  source_eta -= 2*f_r ;
231  resu.set_vr_eta_mu(f_r, source_eta.poisson_angu(), mu().poisson());
232 
233  break ;
234 
235  }
236 
237  case 2 : {
238 
239  Tenseur source_p(*mp, 1, CON, mp->get_bvect_spher() ) ;
240  source_p.set_etat_qcq() ;
241  for (int i=0; i<3; i++) {
242  source_p.set(i) = Cmp(*cmp[i]) ;
243  }
244  source_p.change_triad(mp->get_bvect_cart()) ;
245  Tenseur vect_auxi (*mp, 1, CON, mp->get_bvect_cart()) ;
246  vect_auxi.set_etat_qcq() ;
247  Tenseur scal_auxi (*mp) ;
248  scal_auxi.set_etat_qcq() ;
249 
250  Tenseur resu_p(source_p.poisson_vect(lambda, vect_auxi, scal_auxi)) ;
251  resu_p.change_triad(mp->get_bvect_spher() ) ;
252 
253  for (int i=1; i<=3; i++)
254  resu.set(i) = resu_p(i-1) ;
255 
256  break ;
257  }
258 
259  case 3 : {
260 
261  Tenseur source_p(*mp, 1, CON, mp->get_bvect_spher() ) ;
262  source_p.set_etat_qcq() ;
263  for (int i=0; i<3; i++) {
264  source_p.set(i) = Cmp(*cmp[i]) ;
265  }
266  source_p.change_triad(mp->get_bvect_cart()) ;
267  Tenseur scal_auxi (*mp) ;
268  scal_auxi.set_etat_qcq() ;
269 
270  Tenseur resu_p(source_p.poisson_vect_oohara(lambda, scal_auxi)) ;
271  resu_p.change_triad(mp->get_bvect_spher() ) ;
272 
273  for (int i=1; i<=3; i++)
274  resu.set(i) = resu_p(i-1) ;
275 
276  break ;
277  }
278 
279  case 4 : {
280  Scalar divf(*mp) ;
281  if (fabs(lambda+1) < 1.e-8)
282  divf.set_etat_zero() ;
283  else {
284  divf = (potential(met_f) / (lambda + 1)) ;
285  }
286 
287  int nz = mp->get_mg()->get_nzone() ;
288 
289  //-----------------------------------
290  // Removal of the l=0 part of div(F)
291  //-----------------------------------
292  Scalar div_tmp = divf ;
293  div_tmp.div_r_dzpuis(4) ;
294  Scalar source_r(*mp) ;
295  Valeur& va_div = div_tmp.set_spectral_va() ;
296  if (div_tmp.get_etat() != ETATZERO) {
297  va_div.ylm() ;
298  for (int lz=0; lz<nz; lz++) {
299  int np = mp->get_mg()->get_np(lz) ;
300  int nt = mp->get_mg()->get_nt(lz) ;
301  int nr = mp->get_mg()->get_nr(lz) ;
302  if (va_div.c_cf->operator()(lz).get_etat() != ETATZERO)
303  for (int k=0; k<np+1; k++)
304  for (int j=0; j<nt; j++) {
305  int l_q, m_q, base_r ;
306  if (nullite_plm(j, nt, k, np, va_div.base) == 1) {
307  donne_lm(nz, lz, j, k, va_div.base, m_q, l_q, base_r) ;
308  if (l_q == 0)
309  for (int i=0; i<nr; i++)
310  va_div.c_cf->set(lz, k, j, i) = 0. ;
311  }
312  }
313  }
314  source_r.set_etat_qcq() ;
315  source_r.set_spectral_va() = 2*(*va_div.c_cf) ; //2*div(F)
316  source_r.set_spectral_va().ylm_i() ;
317  source_r.set_dzpuis(4) ;
318  }
319  else
320  source_r.set_etat_zero() ;
321 
322  //------------------------
323  // Other source terms ....
324  //------------------------
325  source_r += *(cmp[0]) - lambda*divf.dsdr() ;
326  Scalar f_r(*mp) ;
327  if (source_r.get_etat() != ETATZERO) {
328 
329  //------------------------
330  // The elliptic operator
331  //------------------------
332  Param_elliptic param_fr(source_r) ;
333  for (int lz=0; lz<nz; lz++)
334  param_fr.set_poisson_vect_r(lz) ;
335 
336  f_r = source_r.sol_elliptic(param_fr) ;
337  }
338  else
339  f_r.set_etat_zero() ;
340 
341  //--------------------------
342  // Equation for eta
343  //--------------------------
344  Scalar source_eta = *cmp[1] ;
345  source_eta.mult_sint() ;
346  source_eta = source_eta.dsdt() ;
347  source_eta.div_sint() ;
348  source_eta = (source_eta + cmp[2]->stdsdp()).poisson_angu() ;
349 
350  Scalar dfrsr = f_r.dsdr() ;
351  dfrsr.div_r_dzpuis(4) ;
352  Scalar frsr2 = f_r ;
353  frsr2.div_r_dzpuis(2) ;
354  frsr2.div_r_dzpuis(4) ;
355 
356  Valeur& va_dfr = dfrsr.set_spectral_va() ;
357  Valeur& va_fsr = frsr2.set_spectral_va() ;
358  va_dfr.ylm() ;
359  va_fsr.ylm() ;
360 
361  //Since the operator depends on the domain, the source
362  //must be transformed accordingly.
363  Valeur& va_eta = source_eta.set_spectral_va() ;
364  if (source_eta.get_etat() == ETATZERO) source_eta.annule_hard() ;
365  va_eta.ylm() ;
366  for (int lz=0; lz<nz; lz++) {
367  bool ced = (mp->get_mg()->get_type_r(lz) == UNSURR) ;
368  int np = mp->get_mg()->get_np(lz) ;
369  int nt = mp->get_mg()->get_nt(lz) ;
370  int nr = mp->get_mg()->get_nr(lz) ;
371  if (va_eta.c_cf->operator()(lz).get_etat() != ETATZERO)
372  for (int k=0; k<np+1; k++)
373  for (int j=0; j<nt; j++) {
374  int l_q, m_q, base_r ;
375  if (nullite_plm(j, nt, k, np, va_eta.base) == 1) {
376  donne_lm(nz, lz, j, k, va_eta.base, m_q, l_q, base_r) ;
377  if (l_q > 0)
378  for (int i=0; i<nr; i++) {
379  if (va_div.c_cf->operator()(lz).get_etat() != ETATZERO)
380  va_eta.c_cf->set(lz, k, j, i)
381  -= (lambda + 2. / double(ced ? -l_q : (l_q+1) ))
382  * va_div.c_cf->operator()(lz, k, j, i) ;
383  if (va_fsr.c_cf->operator()(lz).get_etat() != ETATZERO) {
384  va_eta.c_cf->set(lz, k, j, i)
385  += 2. / double(ced ? -l_q : (l_q+1) )
386  * va_dfr.c_cf->operator()(lz, k, j, i) ;
387  va_eta.c_cf->set(lz, k, j, i)
388  -= 2.*( ced ? double(l_q+2)/double(l_q)
389  : double(l_q-1)/double(l_q+1) )
390  * va_fsr.c_cf->set(lz, k, j, i) ;
391  }
392  } //Loop on r
393  } //nullite_plm
394  } //Loop on theta
395  } // Loop on domains
396  if (va_eta.c != 0x0) {
397  delete va_eta.c;
398  va_eta.c = 0x0 ;
399  }
400  va_eta.ylm_i() ;
401 
402  //------------------------
403  // The elliptic operator
404  //------------------------
405  Param_elliptic param_eta(source_eta) ;
406  for (int lz=0; lz<nz; lz++)
407  param_eta.set_poisson_vect_eta(lz) ;
408 
409  resu.set_vr_eta_mu(f_r, source_eta.sol_elliptic(param_eta), mu().poisson()) ;
410  break ;
411 
412  }
413 
414  case 5 : {
415 
416  Scalar divf(*mp) ;
417  if (fabs(lambda+1) < 1.e-8)
418  divf.set_etat_zero() ;
419  else {
420  divf = (potential(met_f) / (lambda + 1)) ;
421  }
422 
423  int nz = mp->get_mg()->get_nzone() ;
424 
425  //-----------------------------------
426  // Removal of the l=0 part of div(F)
427  //-----------------------------------
428  Scalar div_lnot0 = divf ;
429  div_lnot0.div_r_dzpuis(4) ;
430  Scalar source_r(*mp) ;
431  Valeur& va_div = div_lnot0.set_spectral_va() ;
432  if (div_lnot0.get_etat() != ETATZERO) {
433  va_div.coef() ;
434  va_div.ylm() ;
435  for (int lz=0; lz<nz; lz++) {
436  int np = mp->get_mg()->get_np(lz) ;
437  int nt = mp->get_mg()->get_nt(lz) ;
438  int nr = mp->get_mg()->get_nr(lz) ;
439  if (va_div.c_cf->operator()(lz).get_etat() != ETATZERO)
440  for (int k=0; k<np+1; k++)
441  for (int j=0; j<nt; j++) {
442  int l_q, m_q, base_r ;
443  if (nullite_plm(j, nt, k, np, va_div.base) == 1) {
444  donne_lm(nz, lz, j, k, va_div.base, m_q, l_q, base_r) ;
445  if (l_q == 0)
446  for (int i=0; i<nr; i++)
447  va_div.c_cf->set(lz, k, j, i) = 0. ;
448  }
449  }
450  }
451  source_r.set_etat_qcq() ;
452  source_r.set_spectral_va() = 2*(*va_div.c_cf) ; //2*div(F)
453  source_r.set_spectral_va().ylm_i() ;
454  source_r.set_dzpuis(4) ;
455  }
456  else
457  source_r.set_etat_zero() ;
458 
459  //------------------------
460  // Other source terms ....
461  //------------------------
462  source_r += *(cmp[0]) - lambda*divf.dsdr() ;
463  Scalar f_r(*mp) ;
464  if (source_r.get_etat() != ETATZERO) {
465 
466  //------------------------
467  // The elliptic operator
468  //------------------------
469 
470  Param_elliptic param_fr(source_r) ;
471  for (int lz=0; lz<nz; lz++)
472  param_fr.set_poisson_vect_r(lz) ;
473 
474  f_r = source_r.sol_elliptic(param_fr) ;
475  }
476  else
477  f_r.set_etat_zero() ;
478 
479  Scalar source_eta = - *(cmp[0]) ;
480  source_eta.mult_r_dzpuis(3) ;
481  source_eta -= (lambda+2.)*divf ;
482  source_eta.dec_dzpuis() ;
483  f_r.set_spectral_va().ylm() ;
484  Scalar tmp = 2*f_r + f_r.lapang() ;
485  tmp.div_r_dzpuis(2) ;
486  source_eta += tmp ;
487  tmp = (1.+lambda)*divf ;
488  tmp.mult_r_dzpuis(0) ;
489  tmp += f_r ;
490  source_eta = source_eta.primr() ;
491  f_r.set_spectral_va().ylm_i() ;
492  resu.set_vr_eta_mu(f_r, (tmp+source_eta).poisson_angu(), mu().poisson()) ;
493  break ;
494 
495  }
496 
497  case 6 : {
498 
499  poisson_block(lambda, resu) ;
500  break ;
501 
502  }
503  default : {
504  cout << "Vector::poisson : unexpected type of method !" << endl
505  << " method = " << method << endl ;
506  abort() ;
507  break ;
508  }
509 
510  } // End of switch
511 
512  } // End of non-null case
513 
514  return resu ;
515 
516 }
517 
518 Vector Vector::poisson(double lambda, int method) const {
519 
520  const Base_vect_spher* tspher = dynamic_cast<const Base_vect_spher*>(triad) ;
521  const Base_vect_cart* tcart = dynamic_cast<const Base_vect_cart*>(triad) ;
522 
523  assert ((tspher != 0x0) || (tcart != 0x0)) ;
524  const Metric_flat* met_f = 0x0 ;
525 
526  if (tspher != 0x0) {
527  assert (tcart == 0x0) ;
528  met_f = &(mp->flat_met_spher()) ;
529  }
530 
531  if (tcart != 0x0) {
532  assert (tspher == 0x0) ;
533  met_f = &(mp->flat_met_cart()) ;
534  }
535 
536  return ( poisson(lambda, *met_f, method) );
537 
538 }
539 
540 // Version with parameters
541 // -----------------------
542 
543 Vector Vector::poisson(const double lambda, Param& par, int method) const {
544 
545 
546  for (int i=0; i<3; i++)
547  assert(cmp[i]->check_dzpuis(4)) ;
548 
549  assert ((method==0) || (method==2)) ;
550 
551  Vector resu(*mp, CON, triad) ;
552 
553  switch (method) {
554 
555  case 0 : {
556 
557  Metric_flat met_local(*mp, *triad) ;
558  int nitermax = par.get_int(0) ;
559  int& niter = par.get_int_mod(0) ;
560  double relax = par.get_double(0) ;
561  double precis = par.get_double(1) ;
562  Cmp& ss_phi = par.get_cmp_mod(0) ;
563  Cmp& ss_khi = par.get_cmp_mod(1) ;
564  Cmp& ss_mu = par.get_cmp_mod(2) ;
565 
566  Param par_phi ;
567  par_phi.add_int(nitermax, 0) ;
568  par_phi.add_int_mod(niter, 0) ;
569  par_phi.add_double(relax, 0) ;
570  par_phi.add_double(precis, 1) ;
571  par_phi.add_cmp_mod(ss_phi, 0) ;
572 
573  Scalar poten(*mp) ;
574  if (fabs(lambda+1) < 1.e-8)
575  poten.set_etat_zero() ;
576  else {
577  Scalar tmp = potential(met_local) / (lambda + 1) ;
578  tmp.inc_dzpuis(2) ;
579  tmp.poisson(par_phi, poten) ;
580  }
581 
582  Vector grad = poten.derive_con(met_local) ;
583  grad.dec_dzpuis(2) ;
584 
585  Param par_free ;
586  par_free.add_int(nitermax, 0) ;
587  par_free.add_int_mod(niter, 0) ;
588  par_free.add_double(relax, 0) ;
589  par_free.add_double(precis, 1) ;
590  par_free.add_cmp_mod(ss_khi, 0) ;
591  par_free.add_cmp_mod(ss_mu, 1) ;
592 
593  return div_free(met_local).poisson(par_free) + grad ;
594  break ;
595  }
596 
597 
598 
599  case 2 : {
600 
601  Tenseur source_p(*mp, 1, CON, mp->get_bvect_spher() ) ;
602  source_p.set_etat_qcq() ;
603  for (int i=0; i<3; i++) {
604  source_p.set(i) = Cmp(*cmp[i]) ;
605  }
606  source_p.change_triad(mp->get_bvect_cart()) ;
607 
608  Tenseur vect_auxi (*mp, 1, CON, mp->get_bvect_cart()) ;
609  vect_auxi.set_etat_qcq() ;
610  for (int i=0; i<3 ;i++){
611  vect_auxi.set(i) = 0. ;
612  }
613  Tenseur scal_auxi (*mp) ;
614  scal_auxi.set_etat_qcq() ;
615  scal_auxi.set().annule_hard() ;
616  scal_auxi.set_std_base() ;
617 
618  Tenseur resu_p(*mp, 1, CON, mp->get_bvect_cart() ) ;
619  resu_p.set_etat_qcq() ;
620  source_p.poisson_vect(lambda, par, resu_p, vect_auxi, scal_auxi) ;
621  resu_p.change_triad(mp->get_bvect_spher() ) ;
622 
623  for (int i=1; i<=3; i++)
624  resu.set(i) = resu_p(i-1) ;
625 
626  break ;
627  }
628 
629  default : {
630  cout << "Vector::poisson : unexpected type of method !" << endl
631  << " method = " << method << endl ;
632 
633 
634  abort() ;
635  break ;
636  }
637 
638  }// End of switch
639  return resu ;
640 
641 }
642 
643 
644 
645 
646 }
void set_poisson_vect_eta(int zone)
Sets the operator to be a regular elliptic operator to solve for the component of the vector Poisson...
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:443
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:246
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131
const Scalar & lapang() const
Returns the angular Laplacian of *this , where .
Definition: scalar_deriv.C:461
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:324
const Vector_divfree & div_free(const Metric &) const
Returns the div-free vector in the Helmholtz decomposition.
Definition: vector.C:504
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
Lorene prototypes.
Definition: app_hor.h:64
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:294
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
const Scalar & dsdt() const
Returns of *this .
Definition: scalar_deriv.C:208
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Flat metric for tensor calculation.
Definition: metric.h:261
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: scalar_pde.C:136
Scalar primr(bool null_infty=true) const
Computes the radial primitive which vanishes for .
Definition: scalar_integ.C:101
void mult_sint()
Multiplication by .
void annule_hard()
Sets the Cmp to zero in a hard way.
Definition: cmp.C:338
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tensor.h:303
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:353
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition: scalar_deriv.C:402
Tensor field of valence 1.
Definition: vector.h:188
Scalar poisson_angu(double lambda=0) const
Solves the (generalized) angular Poisson equation with *this as source.
Definition: scalar_pde.C:200
const Metric_flat & flat_met_cart() const
Returns the flat metric associated with the Cartesian coordinates and with components expressed in th...
Definition: map.C:331
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:808
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition: tensor.C:808
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:380
Cmp & get_cmp_mod(int position=0) const
Returns the reference of a modifiable Cmp stored in the list.
Definition: param.C:1049
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:292
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition: tenseur.C:668
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Vector poisson(double lambda, int method=6) const
Solves the vector Poisson equation with *this as a source.
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
void set_vr_eta_mu(const Scalar &vr_i, const Scalar &eta_i, const Scalar &mu_i)
Defines the components through potentials and (see members p_eta and p_mu ), as well as the compon...
Definition: vector_etamu.C:189
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition: tensor.h:315
Scalar sol_elliptic(Param_elliptic &params) const
Resolution of a general elliptic equation, putting zero at infinity.
Definition: scalar_pde.C:234
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:299
Parameter storage.
Definition: param.h:125
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Definition: param.C:430
const Scalar & stdsdp() const
Returns of *this .
Definition: scalar_deriv.C:238
Vector_divfree poisson() const
Computes the solution of a vectorial Poisson equation with *this as a source: .
This class contains the parameters needed to call the general elliptic solver.
Cartesian vectorial bases (triads).
Definition: base_vect.h:201
Spherical orthonormal vectorial bases (triads).
Definition: base_vect.h:308
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:791
const Scalar & potential(const Metric &) const
Returns the potential in the Helmholtz decomposition.
Definition: vector.C:492
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:315
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written: .
Definition: vector_etamu.C:98
void div_sint()
Division by .
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:497
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition: param.C:361
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
Definition: param.C:1004
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:474
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:295
void div_r_dzpuis(int ced_mult_r)
Division by r everywhere but with the output flag dzpuis set to ced_mult_r .
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
Definition: map.C:321
void mult_r_dzpuis(int ced_mult_r)
Multiplication by r everywhere but with the output flag dzpuis set to ced_mult_r .
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:298
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition: param.C:385
void set_poisson_vect_r(int zone, bool only_l_zero=false)
Sets the operator to in all domains, for ; and to in all domains otherwise.