LORENE
sym_tensor_aux.C
1 /*
2  * Methods of class Sym_tensor linked with auxiliary members (eta, mu, W, X...)
3  *
4  * (see file sym_tensor.h for documentation)
5  *
6  */
7 
8 /*
9  * Copyright (c) 2005 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 as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 
32 /*
33  * $Id: sym_tensor_aux.C,v 1.16 2016/12/05 16:18:17 j_novak Exp $
34  * $Log: sym_tensor_aux.C,v $
35  * Revision 1.16 2016/12/05 16:18:17 j_novak
36  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
37  *
38  * Revision 1.15 2014/10/13 08:53:43 j_novak
39  * Lorene classes and functions now belong to the namespace Lorene.
40  *
41  * Revision 1.14 2014/10/06 15:13:19 j_novak
42  * Modified #include directives to use c++ syntax.
43  *
44  * Revision 1.13 2007/11/27 15:49:51 n_vasset
45  * new function compute_tilde_C in class sym_tensor
46  *
47  * Revision 1.12 2006/10/24 14:52:38 j_novak
48  * The Mtbl corresponding to the physical space is destroyed after the
49  * calculation of tilde_B(tt), to get the updated result.
50  *
51  * Revision 1.11 2006/10/24 13:03:19 j_novak
52  * New methods for the solution of the tensor wave equation. Perhaps, first
53  * operational version...
54  *
55  * Revision 1.10 2006/08/31 12:12:43 j_novak
56  * Correction of a small mistake in a phi loop.
57  *
58  * Revision 1.9 2006/06/28 07:48:26 j_novak
59  * Better treatment of some null cases.
60  *
61  * Revision 1.8 2006/06/12 11:40:07 j_novak
62  * Better memory and spectral base handling for Sym_tensor::compute_tilde_B.
63  *
64  * Revision 1.7 2006/06/12 07:42:29 j_novak
65  * Fields A and tilde{B} are defined only for l>1.
66  *
67  * Revision 1.6 2006/06/12 07:27:20 j_novak
68  * New members concerning A and tilde{B}, dealing with the transverse part of the
69  * Sym_tensor.
70  *
71  * Revision 1.5 2005/09/07 16:47:43 j_novak
72  * Removed method Sym_tensor_trans::T_from_det_one
73  * Modified Sym_tensor::set_auxiliary, so that it takes eta/r and mu/r as
74  * arguments.
75  * Modified Sym_tensor_trans::set_hrr_mu.
76  * Added new protected method Sym_tensor_trans::solve_hrr
77  *
78  * Revision 1.4 2005/04/05 15:38:08 j_novak
79  * Changed the formulas for W and X. There was an error before...
80  *
81  * Revision 1.3 2005/04/05 09:22:15 j_novak
82  * Use of the right formula with poisson_angu(2) for the determination of W and
83  * X.
84  *
85  * Revision 1.2 2005/04/04 15:25:24 j_novak
86  * Added new members www, xxx, ttt and the associated methods.
87  *
88  * Revision 1.1 2005/04/01 14:39:57 j_novak
89  * Methods dealing with auxilliary derived members of the symmetric
90  * tensor (eta, mu, W, X, etc ...).
91  *
92  *
93  *
94  * $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_aux.C,v 1.16 2016/12/05 16:18:17 j_novak Exp $
95  *
96  */
97 
98 // Headers C
99 #include <cstdlib>
100 #include <cassert>
101 #include <cmath>
102 
103 // Headers Lorene
104 #include "metric.h"
105 #include "proto.h"
106 
107 
108  //--------------//
109  // eta //
110  //--------------//
111 
112 
113 namespace Lorene {
114 const Scalar& Sym_tensor::eta(Param* par) const {
115 
116  if (p_eta == 0x0) { // a new computation is necessary
117 
118  // All this has a meaning only for spherical components:
119  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
120 
121  // eta is computed from its definition (148) of Bonazzola et al. (2004)
122  int dzp = operator()(1,1).get_dzpuis() ;
123  int dzp_resu = ((dzp == 0) ? 0 : dzp-1) ;
124 
125  Scalar source_eta = operator()(1,2) ;
126  source_eta.div_tant() ;
127  source_eta += operator()(1,2).dsdt() + operator()(1,3).stdsdp() ;
128  source_eta.mult_r_dzpuis(dzp_resu) ;
129 
130  // Resolution of the angular Poisson equation for eta
131  // --------------------------------------------------
132  if (dynamic_cast<const Map_af*>(mp) != 0x0) {
133  p_eta = new Scalar( source_eta.poisson_angu() ) ;
134  }
135  else {
136  Scalar resu (*mp) ;
137  resu = 0. ;
138  mp->poisson_angu(source_eta, *par, resu) ;
139  p_eta = new Scalar( resu ) ;
140  }
141 
142  }
143 
144  return *p_eta ;
145 
146 }
147 
148 
149  //--------------//
150  // mu //
151  //--------------//
152 
153 
154 const Scalar& Sym_tensor::mu(Param* par) const {
155 
156  if (p_mu == 0x0) { // a new computation is necessary
157 
158  // All this has a meaning only for spherical components:
159  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
160 
161  Scalar source_mu = operator()(1,3) ; // h^{r ph}
162  source_mu.div_tant() ; // h^{r ph} / tan(th)
163 
164  // dh^{r ph}/dth + h^{r ph}/tan(th) - 1/sin(th) dh^{r th}/dphi
165  source_mu += operator()(1,3).dsdt() - operator()(1,2).stdsdp() ;
166 
167  // Multiplication by r
168  int dzp = operator()(1,2).get_dzpuis() ;
169  int dzp_resu = ((dzp == 0) ? 0 : dzp-1) ;
170  source_mu.mult_r_dzpuis(dzp_resu) ;
171 
172  // Resolution of the angular Poisson equation for mu
173  // --------------------------------------------------
174  if (dynamic_cast<const Map_af*>(mp) != 0x0) {
175  p_mu = new Scalar( source_mu.poisson_angu() ) ;
176  }
177  else {
178  Scalar resu (*mp) ;
179  resu = 0. ;
180  mp->poisson_angu(source_mu, *par, resu) ;
181  p_mu = new Scalar( resu ) ;
182  }
183  }
184  return *p_mu ;
185 
186 }
187 
188  //-------------//
189  // T //
190  //-------------//
191 
192 
193 const Scalar& Sym_tensor::ttt() const {
194 
195  if (p_ttt == 0x0) { // a new computation is necessary
196 
197  // All this has a meaning only for spherical components:
198  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
199 
200  p_ttt = new Scalar( operator()(2,2) + operator()(3,3) ) ;
201 
202  }
203  return *p_ttt ;
204 
205 }
206 
207  //------------//
208  // W //
209  //------------//
210 
211 
212 const Scalar& Sym_tensor::www() const {
213 
214  if (p_www == 0x0) { // a new computation is necessary
215 
216  // All this has a meaning only for spherical components:
217  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
218 
219  Scalar ppp = 0.5*(operator()(2,2) - operator()(3,3)) ;
220  Scalar tmp = ppp.dsdt() ;
221  tmp.div_tant() ;
222  Scalar source_w = ppp.dsdt().dsdt() +3*tmp - ppp.stdsdp().stdsdp() -2*ppp ;
223  tmp = operator()(2,3) ;
224  tmp.div_tant() ;
225  tmp += operator()(2,3).dsdt() ;
226  source_w += 2*tmp.stdsdp() ;
227 
228  // Resolution of the angular Poisson equation for W
229  p_www = new Scalar( source_w.poisson_angu(2).poisson_angu() ) ;
230 
231  }
232 
233  return *p_www ;
234 
235 }
236 
237 
238  //------------//
239  // X //
240  //------------//
241 
242 
243 const Scalar& Sym_tensor::xxx() const {
244 
245  if (p_xxx == 0x0) { // a new computation is necessary
246 
247  // All this has a meaning only for spherical components:
248  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
249 
250  const Scalar& htp = operator()(2,3) ;
251  Scalar tmp = htp.dsdt() ;
252  tmp.div_tant() ;
253  Scalar source_x = htp.dsdt().dsdt() + 3*tmp - htp.stdsdp().stdsdp() -2*htp ;
254  Scalar ppp = 0.5*(operator()(2,2) - operator()(3,3)) ;
255  tmp = ppp ;
256  tmp.div_tant() ;
257  tmp += ppp.dsdt() ;
258  source_x -= 2*tmp.stdsdp() ;
259 
260  // Resolution of the angular Poisson equation for W
261  p_xxx = new Scalar( source_x.poisson_angu(2).poisson_angu() ) ;
262 
263  }
264 
265  return *p_xxx ;
266 
267 }
268 
269 void Sym_tensor::set_auxiliary(const Scalar& trr, const Scalar& eta_over_r,
270  const Scalar& mu_over_r, const Scalar& w_in,
271  const Scalar& x_in, const Scalar& t_in ) {
272 
273  // All this has a meaning only for spherical components:
274  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
275  int dzp = trr.get_dzpuis() ;
276  int dzeta = (dzp == 0 ? 0 : dzp - 1) ;
277 
278  assert(eta_over_r.check_dzpuis(dzp)) ;
279  assert(mu_over_r.check_dzpuis(dzp)) ;
280  assert(w_in.check_dzpuis(dzp)) ;
281  assert(x_in.check_dzpuis(dzp)) ;
282  assert(t_in.check_dzpuis(dzp)) ;
283 
284  assert(&w_in != p_www) ;
285  assert(&x_in != p_xxx) ;
286  assert(&t_in != p_ttt) ;
287 
288  set(1,1) = trr ;
289  set(1,2) = eta_over_r.dsdt() - mu_over_r.stdsdp() ;
290  // set(1,2).div_r_dzpuis(dzp) ;
291  set(1,3) = eta_over_r.stdsdp() + mu_over_r.dsdt() ;
292  // set(1,3).div_r_dzpuis(dzp) ;
293  Scalar tmp = w_in.lapang() ;
294  tmp.set_spectral_va().ylm_i() ;
295  Scalar ppp = 2*w_in.dsdt().dsdt() - tmp - 2*x_in.stdsdp().dsdt() ;
296  tmp = x_in.lapang() ;
297  tmp.set_spectral_va().ylm_i() ;
298  set(2,3) = 2*x_in.dsdt().dsdt() - tmp + 2*w_in.stdsdp().dsdt() ;
299  set(2,2) = 0.5*t_in + ppp ;
300  set(3,3) = 0.5*t_in - ppp ;
301 
302  // Deleting old derived quantities ...
303  del_deriv() ;
304 
305  // .. and affecting new ones.
306  p_eta = new Scalar(eta_over_r) ; p_eta->mult_r_dzpuis(dzeta) ;
307  p_mu = new Scalar(mu_over_r) ; p_mu->mult_r_dzpuis(dzeta) ;
308  p_www = new Scalar(w_in) ;
309  p_xxx = new Scalar(x_in) ;
310  p_ttt = new Scalar(t_in) ;
311 
312 }
313 
314  //------------//
315  // A //
316  //------------//
317 
318 
319 const Scalar& Sym_tensor::compute_A(bool output_ylm, Param* par) const {
320 
321  if (p_aaa == 0x0) { // a new computation is necessary
322 
323  // All this has a meaning only for spherical components:
324  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
325 
326  int dzp = xxx().get_dzpuis() ;
327  int dzp_resu = ((dzp == 0) ? 2 : dzp+1) ;
328 
329  Scalar source_mu = operator()(1,3) ; // h^{r ph}
330  source_mu.div_tant() ; // h^{r ph} / tan(th)
331 
332  // dh^{r ph}/dth + h^{r ph}/tan(th) - 1/sin(th) dh^{r th}/dphi
333  source_mu += operator()(1,3).dsdt() - operator()(1,2).stdsdp() ;
334 
335  // Resolution of the angular Poisson equation for mu / r
336  // -----------------------------------------------------
337  Scalar tilde_mu(*mp) ;
338  tilde_mu = 0 ;
339 
340  if (dynamic_cast<const Map_af*>(mp) != 0x0) {
341  tilde_mu = source_mu.poisson_angu() ;
342  }
343  else {
344  assert(par != 0x0) ;
345  mp->poisson_angu(source_mu, *par, tilde_mu) ;
346  }
347  tilde_mu.div_r_dzpuis(dzp_resu) ;
348 
349  p_aaa = new Scalar( xxx().dsdr() - tilde_mu) ;
350  p_aaa->annule_l(0, 1, output_ylm) ;
351  }
352 
353  if (output_ylm) p_aaa->set_spectral_va().ylm() ;
354  else p_aaa->set_spectral_va().ylm_i() ;
355 
356  return *p_aaa ;
357 
358 }
359 
360  //--------------//
361  // tilde(B) //
362  //--------------//
363 
364 
365 const Scalar& Sym_tensor::compute_tilde_B(bool output_ylm, Param* par) const {
366 
367  if (p_tilde_b == 0x0) { // a new computation is necessary
368 
369  // All this has a meaning only for spherical components:
370  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
371 
372  int dzp = operator()(1,1).get_dzpuis() ;
373  int dzp_resu = ((dzp == 0) ? 2 : dzp+1) ;
374 
375  p_tilde_b = new Scalar(*mp) ;
377 
378  Scalar source_eta = operator()(1,2) ;
379  source_eta.div_tant() ;
380  source_eta += operator()(1,2).dsdt() + operator()(1,3).stdsdp() ;
381 
382  // Resolution of the angular Poisson equation for eta / r
383  // ------------------------------------------------------
384  Scalar tilde_eta(*mp) ;
385  tilde_eta = 0 ;
386 
387  if (dynamic_cast<const Map_af*>(mp) != 0x0) {
388  tilde_eta = source_eta.poisson_angu() ;
389  }
390  else {
391  assert (par != 0x0) ;
392  mp->poisson_angu(source_eta, *par, tilde_eta) ;
393  }
394 
395  Scalar dwdr = www().dsdr() ;
396  dwdr.set_spectral_va().ylm() ;
397  Scalar wsr = www() ;
398  wsr.div_r_dzpuis(dzp_resu) ;
399  wsr.set_spectral_va().ylm() ;
400  Scalar etasr2 = tilde_eta ;
401  etasr2.div_r_dzpuis(dzp_resu) ;
402  etasr2.set_spectral_va().ylm() ;
403  Scalar dtdr = ttt().dsdr() ;
404  dtdr.set_spectral_va().ylm() ;
405  Scalar tsr = ttt() ;
406  tsr.div_r_dzpuis(dzp_resu) ;
407  tsr.set_spectral_va().ylm() ;
408  Scalar hrrsr = operator()(1,1) ;
409  hrrsr.div_r_dzpuis(dzp_resu) ;
410  hrrsr.set_spectral_va().ylm() ;
411 
412  int nz = mp->get_mg()->get_nzone() ;
413 
414  Base_val base(nz) ;
415 
416  if (wsr.get_etat() != ETATZERO) {
417  base = wsr.get_spectral_base() ;
418  }
419  else {
420  if (etasr2.get_etat() != ETATZERO) {
421  base = etasr2.get_spectral_base() ;
422  }
423  else {
424  if (tsr.get_etat() != ETATZERO) {
425  base = tsr.get_spectral_base() ;
426  }
427  else {
428  if (hrrsr.get_etat() != ETATZERO) {
429  base = hrrsr.get_spectral_base() ;
430  }
431  else { //Everything is zero...
433  return *p_tilde_b ;
434  }
435  }
436  }
437  }
438 
442 
443  if (dwdr.get_etat() == ETATZERO) dwdr.annule_hard() ;
444  if (wsr.get_etat() == ETATZERO) wsr.annule_hard() ;
445  if (etasr2.get_etat() == ETATZERO) etasr2.annule_hard() ;
446  if (dtdr.get_etat() == ETATZERO) dtdr.annule_hard() ;
447  if (tsr.get_etat() == ETATZERO) tsr.annule_hard() ;
448  if (hrrsr.get_etat() == ETATZERO) hrrsr.annule_hard() ;
449 
450  int m_q, l_q, base_r ;
451  for (int lz=0; lz<nz; lz++) {
452  int np = mp->get_mg()->get_np(lz) ;
453  int nt = mp->get_mg()->get_nt(lz) ;
454  int nr = mp->get_mg()->get_nr(lz) ;
455  for (int k=0; k<np+1; k++)
456  for (int j=0; j<nt; j++) {
457  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
458  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1))
459  {
460  for (int i=0; i<nr; i++)
461  p_tilde_b->set_spectral_va().c_cf->set(lz, k, j, i)
462  = (l_q+2)*(*dwdr.get_spectral_va().c_cf)(lz, k, j, i)
463  + l_q*(l_q+2)*(*wsr.get_spectral_va().c_cf)(lz, k, j, i)
464  - 2*(*etasr2.get_spectral_va().c_cf)(lz, k, j, i)
465  + 0.5*double(l_q+2)/double(l_q+1)*(*tsr.get_spectral_va().c_cf)(lz, k, j, i)
466  + 0.5/double(l_q+1)*(*dtdr.get_spectral_va().c_cf)(lz, k, j, i)
467  - 1./double(l_q+1)*(*hrrsr.get_spectral_va().c_cf)(lz, k, j, i) ;
468  }
469  }
470  }
471  p_tilde_b->set_dzpuis(dzp_resu) ;
472  } //End of p_tilde_b == 0x0
473 
474  if (output_ylm) p_tilde_b->set_spectral_va().ylm() ;
475  else p_tilde_b->set_spectral_va().ylm_i() ;
476 
477  return *p_tilde_b ;
478 
479 }
480 
481 Scalar Sym_tensor::compute_tilde_B_tt(bool output_ylm, Param* par) const {
482 
483  Scalar resu = compute_tilde_B(true, par) ;
484  if (resu.get_etat() == ETATZERO) return resu ;
485  else {
486  assert(resu.get_etat() == ETATQCQ) ;
487  assert(resu.get_spectral_va().c_cf != 0x0) ;
488  int dzp = operator()(1,1).get_dzpuis() ;
489  int dzp_resu = ((dzp == 0) ? 2 : dzp+1) ;
490 
491  Scalar hsr = operator()(1,1) + ttt() ;
492  if (hsr.get_etat() == ETATZERO) return resu ;
493  Scalar dhdr = hsr.dsdr() ;
494  dhdr.set_spectral_va().ylm() ;
495  hsr.div_r_dzpuis(dzp_resu) ;
496  hsr.set_spectral_va().ylm() ;
497 
498  int nz = mp->get_mg()->get_nzone() ;
499 
500  const Base_val& base = resu.get_spectral_base() ;
501 
502  if (dhdr.get_etat() == ETATZERO) dhdr.annule_hard() ;
503 
504  int m_q, l_q, base_r ;
505  for (int lz=0; lz<nz; lz++) {
506  int np = mp->get_mg()->get_np(lz) ;
507  int nt = mp->get_mg()->get_nt(lz) ;
508  int nr = mp->get_mg()->get_nr(lz) ;
509  for (int k=0; k<np+1; k++)
510  for (int j=0; j<nt; j++) {
511  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
512  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1))
513  {
514  for (int i=0; i<nr; i++)
515  resu.set_spectral_va().c_cf->set(lz, k, j, i)
516  -= 0.5*(*hsr.get_spectral_va().c_cf)(lz, k, j, i)
517  + 0.5/double(l_q+1)*(*dhdr.get_spectral_va().c_cf)(lz, k, j, i);
518  }
519  }
520  }
521  resu.set_dzpuis(dzp_resu) ;
522  if (resu.set_spectral_va().c != 0x0)
523  delete resu.set_spectral_va().c ;
524  resu.set_spectral_va().c = 0x0 ;
525  } //End of resu != 0
526 
527  if (output_ylm) resu.set_spectral_va().ylm() ;
528  else resu.set_spectral_va().ylm_i() ;
529 
530  return resu ;
531 
532 }
533 
535  tras) const {
536 
537  // All this has a meaning only for spherical components:
538  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
539 
540  Scalar resu = tbtt ;
541  if (resu.get_etat() == ETATZERO) {
542  if (tras.get_etat() == ETATZERO) return resu ;
543  else {
544  resu.annule_hard() ;
545  Base_val base = tras.get_spectral_base() ;
546  base.mult_x() ;
547  resu.set_spectral_base(base) ;
548  }
549  }
550  resu.set_spectral_va().ylm() ;
551  int dzp = operator()(1,1).get_dzpuis() ;
552  int dzp_resu = ((dzp == 0) ? 2 : dzp+1) ;
553 
554  Scalar hsr = tras ;
555  if (hsr.get_etat() == ETATZERO) return resu ;
556  else {
557  Scalar dhdr = hsr.dsdr() ;
558  if (dhdr.get_etat() == ETATZERO) dhdr.annule_hard() ;
559  dhdr.set_spectral_va().ylm() ;
560  hsr.div_r_dzpuis(dzp_resu) ;
561  hsr.set_spectral_va().ylm() ;
562 
563  int nz = mp->get_mg()->get_nzone() ;
564  const Base_val& base = resu.get_spectral_base() ;
565 
566  int m_q, l_q, base_r ;
567  for (int lz=0; lz<nz; lz++) {
568  if ((*resu.get_spectral_va().c_cf)(lz).get_etat() == ETATZERO)
569  resu.get_spectral_va().c_cf->set(lz).annule_hard() ;
570  int np = mp->get_mg()->get_np(lz) ;
571  int nt = mp->get_mg()->get_nt(lz) ;
572  int nr = mp->get_mg()->get_nr(lz) ;
573  for (int k=0; k<np+1; k++)
574  for (int j=0; j<nt; j++) {
575  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
576  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1))
577  {
578  for (int i=0; i<nr; i++)
579  resu.set_spectral_va().c_cf->set(lz, k, j, i)
580  += 0.5*(*hsr.get_spectral_va().c_cf)(lz, k, j, i)
581  + 0.5/double(l_q+1)*(*dhdr.get_spectral_va().c_cf)(lz, k, j, i);
582  }
583  }
584  }
585  resu.set_dzpuis(dzp_resu) ;
586  if (resu.set_spectral_va().c != 0x0)
587  delete resu.set_spectral_va().c ;
588  resu.set_spectral_va().c = 0x0 ;
589  } //End of trace != 0
590  return resu ;
591 }
592 
593 
594  //--------------//
595  // tilde(C) //
596  //--------------//
597 
598 
599 const Scalar& Sym_tensor::compute_tilde_C(bool output_ylm, Param* par) const {
600 
601  if (p_tilde_c == 0x0) { // a new computation is necessary
602 
603  // All this has a meaning only for spherical components:
604  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
605 
606  int dzp = operator()(1,1).get_dzpuis() ;
607  int dzp_resu = ((dzp == 0) ? 2 : dzp+1) ;
608 
609  p_tilde_c = new Scalar(*mp) ;
611 
612  Scalar source_eta = operator()(1,2) ;
613  source_eta.div_tant() ;
614  source_eta += operator()(1,2).dsdt() + operator()(1,3).stdsdp() ;
615 
616  // Resolution of the angular Poisson equation for eta / r
617  // ------------------------------------------------------
618  Scalar tilde_eta(*mp) ;
619  tilde_eta = 0 ;
620 
621  if (dynamic_cast<const Map_af*>(mp) != 0x0) {
622  tilde_eta = source_eta.poisson_angu() ;
623  }
624  else {
625  assert (par != 0x0) ;
626  mp->poisson_angu(source_eta, *par, tilde_eta) ;
627  }
628 
629  Scalar dwdr = www().dsdr() ;
630  dwdr.set_spectral_va().ylm() ;
631  Scalar wsr = www() ;
632  wsr.div_r_dzpuis(dzp_resu) ;
633  wsr.set_spectral_va().ylm() ;
634  Scalar etasr2 = tilde_eta ;
635  etasr2.div_r_dzpuis(dzp_resu) ;
636  etasr2.set_spectral_va().ylm() ;
637  Scalar dtdr = ttt().dsdr() ;
638  dtdr.set_spectral_va().ylm() ;
639  Scalar tsr = ttt() ;
640  tsr.div_r_dzpuis(dzp_resu) ;
641  tsr.set_spectral_va().ylm() ;
642  Scalar hrrsr = operator()(1,1) ;
643  hrrsr.div_r_dzpuis(dzp_resu) ;
644  hrrsr.set_spectral_va().ylm() ;
645 
646  int nz = mp->get_mg()->get_nzone() ;
647 
648  Base_val base(nz) ;
649 
650  if (wsr.get_etat() != ETATZERO) {
651  base = wsr.get_spectral_base() ;
652  }
653  else {
654  if (etasr2.get_etat() != ETATZERO) {
655  base = etasr2.get_spectral_base() ;
656  }
657  else {
658  if (tsr.get_etat() != ETATZERO) {
659  base = tsr.get_spectral_base() ;
660  }
661  else {
662  if (hrrsr.get_etat() != ETATZERO) {
663  base = hrrsr.get_spectral_base() ;
664  }
665  else { //Everything is zero...
667  return *p_tilde_c ;
668  }
669  }
670  }
671  }
672 
676 
677  if (dwdr.get_etat() == ETATZERO) dwdr.annule_hard() ;
678  if (wsr.get_etat() == ETATZERO) wsr.annule_hard() ;
679  if (etasr2.get_etat() == ETATZERO) etasr2.annule_hard() ;
680  if (dtdr.get_etat() == ETATZERO) dtdr.annule_hard() ;
681  if (tsr.get_etat() == ETATZERO) tsr.annule_hard() ;
682  if (hrrsr.get_etat() == ETATZERO) hrrsr.annule_hard() ;
683 
684  int m_q, l_q, base_r ;
685  for (int lz=0; lz<nz; lz++) {
686  int np = mp->get_mg()->get_np(lz) ;
687  int nt = mp->get_mg()->get_nt(lz) ;
688  int nr = mp->get_mg()->get_nr(lz) ;
689  for (int k=0; k<np+1; k++)
690  for (int j=0; j<nt; j++) {
691  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
692  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1))
693  {
694  for (int i=0; i<nr; i++)
695  p_tilde_c->set_spectral_va().c_cf->set(lz, k, j, i)
696  = - (l_q - 1)*(*dwdr.get_spectral_va().c_cf)(lz, k, j, i)
697  + (l_q + 1)*(l_q - 1)*(*wsr.get_spectral_va().c_cf)(lz, k, j, i)
698  - 2*(*etasr2.get_spectral_va().c_cf)(lz, k, j, i)
699  + 0.5*double(l_q-1)/double(l_q)*(*tsr.get_spectral_va().c_cf)(lz, k, j, i)
700  - 0.5/double(l_q)*(*dtdr.get_spectral_va().c_cf)(lz, k, j, i)
701  + 1./double(l_q)*(*hrrsr.get_spectral_va().c_cf)(lz, k, j, i) ;
702  }
703  }
704  }
705  p_tilde_c->set_dzpuis(dzp_resu) ;
706  } //End of p_tilde_c == 0x0
707 
708  if (output_ylm) p_tilde_c->set_spectral_va().ylm() ;
709  else p_tilde_c->set_spectral_va().ylm_i() ;
710 
711  return *p_tilde_c ;
712 
713 }
714 }
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Definition: scalar.h:1328
Scalar * p_mu
Field such that the components of the tensor are written (has only meaning with spherical component...
Definition: sym_tensor.h:277
void annule_l(int l_min, int l_max, bool ylm_output=false)
Sets all the multipolar components between l_min and l_max to zero.
Definition: scalar_manip.C:117
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
const Scalar & mu(Param *par=0x0) const
Gives the field (see member p_mu ).
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 set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
Definition: valeur.C:715
Scalar * p_tilde_b
Field defined from and h insensitive to the longitudinal part of the Sym_tensor.
Definition: sym_tensor.h:337
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
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
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
const Scalar & xxx() const
Gives the field X (see member p_xxx ).
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tensor.h:309
Scalar * p_www
Field W such that the components and of the tensor are written (has only meaning with spherical com...
Definition: sym_tensor.h:296
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
const Scalar & compute_tilde_C(bool output_ylm=true, Param *par=0x0) const
Gives the field (see member p_tilde_c ).
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:359
Scalar poisson_angu(double lambda=0) const
Solves the (generalized) angular Poisson equation with *this as source.
Definition: scalar_pde.C:203
Scalar get_tilde_B_from_TT_trace(const Scalar &tilde_B_tt_in, const Scalar &trace) const
Computes (see Sym_tensor::p_tilde_b ) from its transverse-traceless part and the trace...
const Scalar & ttt() const
Gives the field T (see member p_ttt ).
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:814
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:386
Scalar * p_eta
Field such that the components of the tensor are written (has only meaning with spherical component...
Definition: sym_tensor.h:263
Scalar * p_aaa
Field A defined from X and insensitive to the longitudinal part of the Sym_tensor (only for )...
Definition: sym_tensor.h:325
const Scalar & www() const
Gives the field W (see member p_www ).
virtual void del_deriv() const
Deletes the derived quantities.
Definition: sym_tensor.C:289
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:563
const Scalar & compute_A(bool output_ylm=true, Param *par=0x0) const
Gives the field A (see member p_aaa ).
Scalar * p_xxx
Field X such that the components and of the tensor are written (has only meaning with spherical com...
Definition: sym_tensor.h:315
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:309
const Scalar & compute_tilde_B(bool output_ylm=true, Param *par=0x0) const
Gives the field (see member p_tilde_b ).
Parameter storage.
Definition: param.h:125
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Scalar * p_tilde_c
Field defined from and h insensitive to the longitudinal part of the Sym_tensor.
Definition: sym_tensor.h:349
const Scalar & stdsdp() const
Returns of *this .
Definition: scalar_deriv.C:238
virtual const Scalar & eta(Param *par=0x0) const
Gives the field (see member p_eta ).
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
Definition: scalar.C:803
void set_auxiliary(const Scalar &trr, const Scalar &eta_over_r, const Scalar &mu_over_r, const Scalar &www, const Scalar &xxx, const Scalar &ttt)
Assigns the component and the derived members p_eta , p_mu , p_www, p_xxx and p_ttt ...
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Scalar * p_ttt
Field T defined as .
Definition: sym_tensor.h:318
Bases of the spectral expansions.
Definition: base_val.h:325
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition: mtbl_cf.C:315
Scalar compute_tilde_B_tt(bool output_ylm=true, Param *par=0x0) const
Gives the field (see member p_tilde_b ) associated with the TT-part of the Sym_tensor ...
void mult_x()
The basis is transformed as with a multiplication by .
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
const Scalar & operator()(const Itbl &ind) const
Returns the value of a component (read-only version).
Definition: tensor.C:807
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
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 .
void div_tant()
Division by .
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:375
virtual void poisson_angu(const Scalar &source, Param &par, Scalar &uu, double lambda=0) const =0
Computes the solution of the generalized angular Poisson equation.
void mult_r_dzpuis(int ced_mult_r)
Multiplication by r everywhere but with the output flag dzpuis set to ced_mult_r ...
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607