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