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