LORENE
sym_tensor_trans_dirac.C
1 /*
2  * Methods to impose the Dirac gauge: divergence-free condition.
3  *
4  * (see file sym_tensor.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2006 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 sym_tensor_trans_dirac_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_trans_dirac.C,v 1.8 2015/08/10 15:32:26 j_novak Exp $" ;
29 
30 /*
31  * $Id: sym_tensor_trans_dirac.C,v 1.8 2015/08/10 15:32:26 j_novak Exp $
32  * $Log: sym_tensor_trans_dirac.C,v $
33  * Revision 1.8 2015/08/10 15:32:26 j_novak
34  * Better calls to Param::add_int(), to avoid weird problems (e.g. with g++ 4.8).
35  *
36  * Revision 1.7 2014/10/13 08:53:43 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.6 2014/10/06 15:13:19 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.5 2008/08/29 13:15:22 j_novak
43  * Corrected a mistake in the case of no CED.
44  *
45  * Revision 1.4 2008/08/29 05:33:18 j_novak
46  * Minor modif.
47  *
48  * Revision 1.3 2008/08/27 08:13:20 j_novak
49  * Correction of a mistake in the imposition of BCs in sol_Dirac_A. + Treatment of
50  * the case of BCs imposed on a nucleus (nz_bc = 0).
51  *
52  * Revision 1.2 2006/10/24 13:03:19 j_novak
53  * New methods for the solution of the tensor wave equation. Perhaps, first
54  * operational version...
55  *
56  * Revision 1.1 2006/09/05 15:38:45 j_novak
57  * The fuctions sol_Dirac... are in a seperate file, with new parameters to
58  * control the boundary conditions.
59  *
60  *
61  * $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_trans_dirac.C,v 1.8 2015/08/10 15:32:26 j_novak Exp $
62  *
63  */
64 
65 // C headers
66 #include <cassert>
67 #include <cmath>
68 
69 // Lorene headers
70 #include "tensor.h"
71 #include "diff.h"
72 #include "proto.h"
73 #include "param.h"
74 
75 //----------------------------------------------------------------------------------
76 //
77 // sol_Dirac_A
78 //
79 //----------------------------------------------------------------------------------
80 
81 namespace Lorene {
82 void Sym_tensor_trans::sol_Dirac_A(const Scalar& aaa, Scalar& tilde_mu, Scalar& x_new,
83  const Param* par_bc) const {
84 
85  const Map_af* mp_aff = dynamic_cast<const Map_af*>(mp) ;
86  assert(mp_aff != 0x0) ; //Only affine mapping for the moment
87 
88  const Mg3d& mgrid = *mp_aff->get_mg() ;
89  assert(mgrid.get_type_r(0) == RARE) ;
90  if (aaa.get_etat() == ETATZERO) {
91  tilde_mu = 0 ;
92  x_new = 0 ;
93  return ;
94  }
95  assert(aaa.get_etat() != ETATNONDEF) ;
96  int nz = mgrid.get_nzone() ;
97  int nzm1 = nz - 1 ;
98  bool ced = (mgrid.get_type_r(nzm1) == UNSURR) ;
99  int n_shell = ced ? nz-2 : nzm1 ;
100  int nz_bc = nzm1 ;
101  if (par_bc != 0x0)
102  if (par_bc->get_n_int() > 0) nz_bc = par_bc->get_int() ;
103  n_shell = (nz_bc < n_shell ? nz_bc : n_shell) ;
104  bool cedbc = (ced && (nz_bc == nzm1)) ;
105 #ifndef NDEBUG
106  if (!cedbc) {
107  assert(par_bc != 0x0) ;
108  assert(par_bc->get_n_tbl_mod() >= 3) ;
109  }
110 #endif
111  int nt = mgrid.get_nt(0) ;
112  int np = mgrid.get_np(0) ;
113 
114  Scalar source = aaa ;
115  Scalar source_coq = aaa ;
116  source_coq.annule_domain(0) ;
117  if (ced) source_coq.annule_domain(nzm1) ;
118  source_coq.mult_r() ;
119  source.set_spectral_va().ylm() ;
120  source_coq.set_spectral_va().ylm() ;
121  Base_val base = source.get_spectral_base() ;
122  base.mult_x() ;
123 
124  tilde_mu.annule_hard() ;
125  tilde_mu.set_spectral_base(base) ;
126  tilde_mu.set_spectral_va().set_etat_cf_qcq() ;
127  tilde_mu.set_spectral_va().c_cf->annule_hard() ;
128  x_new.annule_hard() ;
129  x_new.set_spectral_base(base) ;
130  x_new.set_spectral_va().set_etat_cf_qcq() ;
131  x_new.set_spectral_va().c_cf->annule_hard() ;
132 
133  Mtbl_cf sol_part_mu(mgrid, base) ; sol_part_mu.annule_hard() ;
134  Mtbl_cf sol_part_x(mgrid, base) ; sol_part_x.annule_hard() ;
135  Mtbl_cf sol_hom1_mu(mgrid, base) ; sol_hom1_mu.annule_hard() ;
136  Mtbl_cf sol_hom1_x(mgrid, base) ; sol_hom1_x.annule_hard() ;
137  Mtbl_cf sol_hom2_mu(mgrid, base) ; sol_hom2_mu.annule_hard() ;
138  Mtbl_cf sol_hom2_x(mgrid, base) ; sol_hom2_x.annule_hard() ;
139 
140  int l_q, m_q, base_r ;
141 
142  //---------------
143  //-- NUCLEUS ---
144  //---------------
145  {int lz = 0 ;
146  int nr = mgrid.get_nr(lz) ;
147  double alpha = mp_aff->get_alpha()[lz] ;
148  Matrice ope(2*nr, 2*nr) ;
149  ope.set_etat_qcq() ;
150 
151  for (int k=0 ; k<np+1 ; k++) {
152  for (int j=0 ; j<nt ; j++) {
153  // quantic numbers and spectral bases
154  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
155  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
156  Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
157  Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
158 
159  for (int lin=0; lin<nr; lin++)
160  for (int col=0; col<nr; col++)
161  ope.set(lin,col) = md(lin,col) + 3*ms(lin,col) ;
162  for (int lin=0; lin<nr; lin++)
163  for (int col=0; col<nr; col++)
164  ope.set(lin,col+nr) = (2-l_q*(l_q+1))*ms(lin,col) ;
165  for (int lin=0; lin<nr; lin++)
166  for (int col=0; col<nr; col++)
167  ope.set(lin+nr,col) = -ms(lin,col) ;
168  for (int lin=0; lin<nr; lin++)
169  for (int col=0; col<nr; col++)
170  ope.set(lin+nr,col+nr) = md(lin,col) ;
171 
172  ope *= 1./alpha ;
173  int ind1 = nr ;
174  for (int col=0; col<2*nr; col++)
175  ope.set(ind1+nr-2, col) = 0 ;
176  for (int col=0; col<2*nr; col++) {
177  ope.set(nr-1, col) = 0 ;
178  ope.set(2*nr-1, col) = 0 ;
179  }
180  int pari = 1 ;
181  if (base_r == R_CHEBP) {
182  for (int col=0; col<nr; col++) {
183  ope.set(nr-1, col) = pari ;
184  ope.set(2*nr-1, col+nr) = pari ;
185  pari = - pari ;
186  }
187  }
188  else { //In the odd case, the last coefficient must be zero!
189  ope.set(nr-1, nr-1) = 1 ;
190  ope.set(2*nr-1, 2*nr-1) = 1 ;
191  }
192  ope.set(ind1+nr-2, ind1) = 1 ;
193  ope.set_lu() ;
194 
195  Tbl sec(2*nr) ;
196  sec.set_etat_qcq() ;
197  for (int lin=0; lin<nr; lin++)
198  sec.set(lin) = 0 ;
199  for (int lin=0; lin<nr; lin++)
200  sec.set(nr+lin) = (*source.get_spectral_va().c_cf)
201  (lz, k, j, lin) ;
202  sec.set(2*nr-1) = 0 ;
203  sec.set(ind1+nr-2) = 0 ;
204  Tbl sol = ope.inverse(sec) ;
205  for (int i=0; i<nr; i++) {
206  sol_part_mu.set(lz, k, j, i) = sol(i) ;
207  sol_part_x.set(lz, k, j, i) = sol(i+nr) ;
208  }
209  sec.annule_hard() ;
210  sec.set(ind1+nr-2) = 1 ;
211  sol = ope.inverse(sec) ;
212  for (int i=0; i<nr; i++) {
213  sol_hom2_mu.set(lz, k, j, i) = sol(i) ;
214  sol_hom2_x.set(lz, k, j, i) = sol(i+nr) ;
215  }
216  }
217  }
218  }
219  }
220 
221  //-------------
222  // -- Shells --
223  //-------------
224 
225  for (int lz=1; lz <= n_shell; lz++) {
226  int nr = mgrid.get_nr(lz) ;
227  assert(mgrid.get_nt(lz) == nt) ;
228  assert(mgrid.get_np(lz) == np) ;
229  double alpha = mp_aff->get_alpha()[lz] ;
230  double ech = mp_aff->get_beta()[lz] / alpha ;
231  Matrice ope(2*nr, 2*nr) ;
232  ope.set_etat_qcq() ;
233 
234  for (int k=0 ; k<np+1 ; k++) {
235  for (int j=0 ; j<nt ; j++) {
236  // quantic numbers and spectral bases
237  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
238  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
239  Diff_xdsdx oxd(base_r, nr) ; const Matrice& mxd = oxd.get_matrice() ;
240  Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
241  Diff_id oid(base_r, nr) ; const Matrice& mid = oid.get_matrice() ;
242 
243  for (int lin=0; lin<nr; lin++)
244  for (int col=0; col<nr; col++)
245  ope.set(lin,col) = mxd(lin,col) + ech*md(lin,col)
246  + 3*mid(lin,col) ;
247  for (int lin=0; lin<nr; lin++)
248  for (int col=0; col<nr; col++)
249  ope.set(lin,col+nr) = (2-l_q*(l_q+1))*mid(lin,col) ;
250  for (int lin=0; lin<nr; lin++)
251  for (int col=0; col<nr; col++)
252  ope.set(lin+nr,col) = -mid(lin,col) ;
253  for (int lin=0; lin<nr; lin++)
254  for (int col=0; col<nr; col++)
255  ope.set(lin+nr,col+nr) = mxd(lin,col) + ech*md(lin,col) ;
256 
257  int ind0 = 0 ;
258  int ind1 = nr ;
259  for (int col=0; col<2*nr; col++) {
260  ope.set(ind0+nr-1, col) = 0 ;
261  ope.set(ind1+nr-1, col) = 0 ;
262  }
263  ope.set(ind0+nr-1, ind0) = 1 ;
264  ope.set(ind1+nr-1, ind1) = 1 ;
265 
266  ope.set_lu() ;
267 
268  Tbl sec(2*nr) ;
269  sec.set_etat_qcq() ;
270  for (int lin=0; lin<nr; lin++)
271  sec.set(lin) = 0 ;
272  for (int lin=0; lin<nr; lin++)
273  sec.set(nr+lin) = (*source_coq.get_spectral_va().c_cf)
274  (lz, k, j, lin) ;
275  sec.set(ind0+nr-1) = 0 ;
276  sec.set(ind1+nr-1) = 0 ;
277  Tbl sol = ope.inverse(sec) ;
278  for (int i=0; i<nr; i++) {
279  sol_part_mu.set(lz, k, j, i) = sol(i) ;
280  sol_part_x.set(lz, k, j, i) = sol(i+nr) ;
281  }
282  sec.annule_hard() ;
283  sec.set(ind0+nr-1) = 1 ;
284  sol = ope.inverse(sec) ;
285  for (int i=0; i<nr; i++) {
286  sol_hom1_mu.set(lz, k, j, i) = sol(i) ;
287  sol_hom1_x.set(lz, k, j, i) = sol(i+nr) ;
288  }
289  sec.set(ind0+nr-1) = 0 ;
290  sec.set(ind1+nr-1) = 1 ;
291  sol = ope.inverse(sec) ;
292  for (int i=0; i<nr; i++) {
293  sol_hom2_mu.set(lz, k, j, i) = sol(i) ;
294  sol_hom2_x.set(lz, k, j, i) = sol(i+nr) ;
295  }
296  }
297  }
298  }
299  }
300 
301  //------------------------------
302  // Compactified external domain
303  //------------------------------
304  if (cedbc) {int lz = nzm1 ;
305  int nr = mgrid.get_nr(lz) ;
306  assert(mgrid.get_nt(lz) == nt) ;
307  assert(mgrid.get_np(lz) == np) ;
308  double alpha = mp_aff->get_alpha()[lz] ;
309  Matrice ope(2*nr, 2*nr) ;
310  ope.set_etat_qcq() ;
311 
312  for (int k=0 ; k<np+1 ; k++) {
313  for (int j=0 ; j<nt ; j++) {
314  // quantic numbers and spectral bases
315  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
316  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
317  Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
318  Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
319 
320  for (int lin=0; lin<nr; lin++)
321  for (int col=0; col<nr; col++)
322  ope.set(lin,col) = - md(lin,col) + 3*ms(lin,col) ;
323  for (int lin=0; lin<nr; lin++)
324  for (int col=0; col<nr; col++)
325  ope.set(lin,col+nr) = (2-l_q*(l_q+1))*ms(lin,col) ;
326  for (int lin=0; lin<nr; lin++)
327  for (int col=0; col<nr; col++)
328  ope.set(lin+nr,col) = -ms(lin,col) ;
329  for (int lin=0; lin<nr; lin++)
330  for (int col=0; col<nr; col++)
331  ope.set(lin+nr,col+nr) = -md(lin,col) ;
332 
333  ope *= 1./alpha ;
334  int ind0 = 0 ;
335  int ind1 = nr ;
336  for (int col=0; col<2*nr; col++) {
337  ope.set(ind0+nr-1, col) = 0 ;
338  ope.set(ind1+nr-2, col) = 0 ;
339  ope.set(ind1+nr-1, col) = 0 ;
340  }
341  for (int col=0; col<nr; col++) {
342  ope.set(ind0+nr-1, col+ind0) = 1 ;
343  ope.set(ind1+nr-1, col+ind1) = 1 ;
344  }
345  ope.set(ind1+nr-2, ind1+1) = 1 ;
346 
347  ope.set_lu() ;
348 
349  Tbl sec(2*nr) ;
350  sec.set_etat_qcq() ;
351  for (int lin=0; lin<nr; lin++)
352  sec.set(lin) = 0 ;
353  for (int lin=0; lin<nr; lin++)
354  sec.set(nr+lin) = (*source.get_spectral_va().c_cf)
355  (lz, k, j, lin) ;
356  sec.set(ind0+nr-1) = 0 ;
357  sec.set(ind1+nr-2) = 0 ;
358  sec.set(ind1+nr-1) = 0 ;
359  Tbl sol = ope.inverse(sec) ;
360  for (int i=0; i<nr; i++) {
361  sol_part_mu.set(lz, k, j, i) = sol(i) ;
362  sol_part_x.set(lz, k, j, i) = sol(i+nr) ;
363  }
364  sec.annule_hard() ;
365  sec.set(ind1+nr-2) = 1 ;
366  sol = ope.inverse(sec) ;
367  for (int i=0; i<nr; i++) {
368  sol_hom1_mu.set(lz, k, j, i) = sol(i) ;
369  sol_hom1_x.set(lz, k, j, i) = sol(i+nr) ;
370  }
371  }
372  }
373  }
374  }
375 
376  //---------------------------------------------------------------
377  // Matching of the solutions across the domains and imposition of
378  // boundary conditions (if no compactified domain)
379  //---------------------------------------------------------------
380  int taille = 2*nz_bc + 1;
381  if (cedbc) taille-- ;
382  Mtbl_cf& mmu = *tilde_mu.set_spectral_va().c_cf ;
383  Mtbl_cf& mw = *x_new.set_spectral_va().c_cf ;
384 
385  Tbl sec_membre(taille) ;
386  Matrice systeme(taille, taille) ;
387  int ligne ; int colonne ;
388  Tbl pipo(1) ;
389  const Tbl& mub = (cedbc ? pipo : par_bc->get_tbl_mod(2) );
390  double c_mu = (cedbc ? 0 : par_bc->get_tbl_mod(0)(0) ) ;
391  double d_mu = (cedbc ? 0 : par_bc->get_tbl_mod(0)(1) ) ;
392  double c_x = (cedbc ? 0 : par_bc->get_tbl_mod(0)(2) ) ;
393  double d_x = (cedbc ? 0 : par_bc->get_tbl_mod(0)(3) ) ;
394  Mtbl_cf dhom1_mu = sol_hom1_mu ;
395  Mtbl_cf dhom2_mu = sol_hom2_mu ;
396  Mtbl_cf dpart_mu = sol_part_mu ;
397  Mtbl_cf dhom1_x = sol_hom1_x ;
398  Mtbl_cf dhom2_x = sol_hom2_x ;
399  Mtbl_cf dpart_x = sol_part_x ;
400  if (!cedbc) {
401  dhom1_mu.dsdx() ;
402  dhom2_mu.dsdx() ;
403  dpart_mu.dsdx() ;
404  dhom1_x.dsdx() ;
405  dhom2_x.dsdx() ;
406  dpart_x.dsdx() ;
407  }
408 
409  // Loop on l and m
410  //----------------
411  for (int k=0 ; k<np+1 ; k++)
412  for (int j=0 ; j<nt ; j++) {
413  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
414  if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
415  ligne = 0 ;
416  colonne = 0 ;
417  systeme.annule_hard() ;
418  sec_membre.annule_hard() ;
419 
420  //Nucleus
421  int nr = mgrid.get_nr(0) ;
422 
423  if (nz_bc >0) {
424  systeme.set(ligne, colonne) = sol_hom2_mu.val_out_bound_jk(0, j, k) ;
425  sec_membre.set(ligne) = -sol_part_mu.val_out_bound_jk(0, j, k) ;
426  ligne++ ;
427 
428  systeme.set(ligne, colonne) = sol_hom2_x.val_out_bound_jk(0, j, k) ;
429  sec_membre.set(ligne) = -sol_part_x.val_out_bound_jk(0, j, k) ;
430  colonne++ ;
431  }
432  //shells
433  for (int zone=1 ; zone<nz_bc ; zone++) {
434  nr = mgrid.get_nr(zone) ;
435  ligne-- ;
436 
437  //Condition at x = -1
438  systeme.set(ligne, colonne) =
439  - sol_hom1_mu.val_in_bound_jk(zone, j, k) ;
440  systeme.set(ligne, colonne+1) =
441  - sol_hom2_mu.val_in_bound_jk(zone, j, k) ;
442 
443  sec_membre.set(ligne) += sol_part_mu.val_in_bound_jk(zone, j, k) ;
444  ligne++ ;
445 
446  systeme.set(ligne, colonne) =
447  - sol_hom1_x.val_in_bound_jk(zone, j, k) ;
448  systeme.set(ligne, colonne+1) =
449  - sol_hom2_x.val_in_bound_jk(zone, j, k) ;
450 
451  sec_membre.set(ligne) += sol_part_x.val_in_bound_jk(zone, j, k) ;
452  ligne++ ;
453 
454  // Condition at x=1
455  systeme.set(ligne, colonne) =
456  sol_hom1_mu.val_out_bound_jk(zone, j, k) ;
457  systeme.set(ligne, colonne+1) =
458  sol_hom2_mu.val_out_bound_jk(zone, j, k) ;
459 
460  sec_membre.set(ligne) -= sol_part_mu.val_out_bound_jk(zone, j, k) ;
461  ligne++ ;
462 
463  systeme.set(ligne, colonne) =
464  sol_hom1_x.val_out_bound_jk(zone, j, k) ;
465  systeme.set(ligne, colonne+1) =
466  sol_hom2_x.val_out_bound_jk(zone, j, k) ;
467 
468  sec_membre.set(ligne) -= sol_part_x.val_out_bound_jk(zone, j, k) ;
469 
470  colonne += 2 ;
471  }
472 
473  //Last domain
474  nr = mgrid.get_nr(nz_bc) ;
475  double alpha = mp_aff->get_alpha()[nz_bc] ;
476  if (nz_bc>0) {
477  ligne-- ;
478 
479  //Condition at x = -1
480  systeme.set(ligne, colonne) =
481  - sol_hom1_mu.val_in_bound_jk(nz_bc, j, k) ;
482  if (!cedbc) systeme.set(ligne, colonne+1) =
483  - sol_hom2_mu.val_in_bound_jk(nz_bc, j, k) ;
484 
485  sec_membre.set(ligne) += sol_part_mu.val_in_bound_jk(nz_bc, j, k) ;
486  ligne++ ;
487 
488  systeme.set(ligne, colonne) =
489  - sol_hom1_x.val_in_bound_jk(nz_bc, j, k) ;
490  if (!cedbc) systeme.set(ligne, colonne+1) =
491  - sol_hom2_x.val_in_bound_jk(nz_bc, j, k) ;
492 
493  sec_membre.set(ligne) += sol_part_x.val_in_bound_jk(nz_bc, j, k) ;
494  ligne++ ;
495  }
496  if (!cedbc) {// Special condition at x=1
497  if (nz_bc>0) {
498  systeme.set(ligne, colonne) =
499  c_mu*sol_hom1_mu.val_out_bound_jk(nz_bc, j, k)
500  + d_mu*dhom1_mu.val_out_bound_jk(nz_bc, j, k) / alpha
501  + c_x*sol_hom1_x.val_out_bound_jk(nz_bc, j, k)
502  + d_x*dhom1_x.val_out_bound_jk(nz_bc, j, k) / alpha ;
503  }
504  else {
505  assert(ligne == 0) ;
506  colonne = -1 ;
507  }
508  systeme.set(ligne, colonne+1) =
509  c_mu*sol_hom2_mu.val_out_bound_jk(nz_bc, j, k)
510  + d_mu*dhom2_mu.val_out_bound_jk(nz_bc, j, k) / alpha
511  + c_x*sol_hom2_x.val_out_bound_jk(nz_bc, j, k)
512  + d_x*dhom2_x.val_out_bound_jk(nz_bc, j, k) / alpha ;
513 
514  sec_membre.set(ligne) -= c_mu*sol_part_mu.val_out_bound_jk(nz_bc, j, k)
515  + d_mu*dpart_mu.val_out_bound_jk(nz_bc, j, k)/alpha
516  + c_x*sol_part_x.val_out_bound_jk(nz_bc, j, k)
517  + d_x*dpart_x.val_out_bound_jk(nz_bc, j, k)/alpha
518  - mub(k, j) ;
519  }
520 
521  // Solution of the system giving the coefficients for the homogeneous
522  // solutions
523  //-------------------------------------------------------------------
524  systeme.set_lu() ;
525  Tbl facteur = systeme.inverse(sec_membre) ;
526  int conte = 0 ;
527 
528  // everything is put to the right place...
529  //----------------------------------------
530  nr = mgrid.get_nr(0) ; //nucleus
531  for (int i=0 ; i<nr ; i++) {
532  mmu.set(0, k, j, i) = sol_part_mu(0, k, j, i)
533  + facteur(conte)*sol_hom2_mu(0, k, j, i) ;
534  mw.set(0, k, j, i) = sol_part_x(0, k, j, i)
535  + facteur(conte)*sol_hom2_x(0, k, j, i) ;
536  }
537  conte++ ;
538  for (int zone=1 ; zone<=n_shell ; zone++) { //shells
539  nr = mgrid.get_nr(zone) ;
540  for (int i=0 ; i<nr ; i++) {
541  mmu.set(zone, k, j, i) = sol_part_mu(zone, k, j, i)
542  + facteur(conte)*sol_hom1_mu(zone, k, j, i)
543  + facteur(conte+1)*sol_hom2_mu(zone, k, j, i) ;
544 
545  mw.set(zone, k, j, i) = sol_part_x(zone, k, j, i)
546  + facteur(conte)*sol_hom1_x(zone, k, j, i)
547  + facteur(conte+1)*sol_hom2_x(zone, k, j, i) ;
548  }
549  conte+=2 ;
550  }
551  if (cedbc) {
552  nr = mgrid.get_nr(nzm1) ; //compactified external domain
553  for (int i=0 ; i<nr ; i++) {
554  mmu.set(nzm1, k, j, i) = sol_part_mu(nzm1, k, j, i)
555  + facteur(conte)*sol_hom1_mu(nzm1, k, j, i) ;
556 
557  mw.set(nzm1, k, j, i) = sol_part_x(nzm1, k, j, i)
558  + facteur(conte)*sol_hom1_x(nzm1, k, j, i) ;
559  }
560  }
561 
562  } // End of nullite_plm
563  } //End of loop on theta
564 
565  if (tilde_mu.set_spectral_va().c != 0x0)
566  delete tilde_mu.set_spectral_va().c ;
567  tilde_mu.set_spectral_va().c = 0x0 ;
568  tilde_mu.set_spectral_va().ylm_i() ;
569 
570  if (x_new.set_spectral_va().c != 0x0)
571  delete x_new.set_spectral_va().c ;
572  x_new.set_spectral_va().c = 0x0 ;
573  x_new.set_spectral_va().ylm_i() ;
574 
575 }
576 
577 //----------------------------------------------------------------------------------
578 //
579 // sol_Dirac_tilde_B
580 //
581 //----------------------------------------------------------------------------------
582 
583 void Sym_tensor_trans::sol_Dirac_tilde_B(const Scalar& tilde_b, const Scalar& hh,
584  Scalar& hrr, Scalar& tilde_eta, Scalar& ww,
585  Param* par_bc, Param* par_mat) const {
586 
587  const Map_af* mp_aff = dynamic_cast<const Map_af*>(mp) ;
588  assert(mp_aff != 0x0) ; //Only affine mapping for the moment
589 
590  const Mg3d& mgrid = *mp_aff->get_mg() ;
591  assert(mgrid.get_type_r(0) == RARE) ;
592  if ( (tilde_b.get_etat() == ETATZERO) && (hh.get_etat() == ETATZERO) ) {
593  hrr = 0 ;
594  tilde_eta = 0 ;
595  ww = 0 ;
596  return ;
597  }
598  int nz = mgrid.get_nzone() ;
599  int nzm1 = nz - 1 ;
600  bool ced = (mgrid.get_type_r(nzm1) == UNSURR) ;
601  int n_shell = ced ? nz-2 : nzm1 ;
602  int nz_bc = nzm1 ;
603  if (par_bc != 0x0)
604  if (par_bc->get_n_int() > 0)
605  nz_bc = par_bc->get_int() ;
606  n_shell = (nz_bc < n_shell ? nz_bc : n_shell) ;
607  bool cedbc = (ced && (nz_bc == nzm1)) ;
608 #ifndef NDEBUG
609  if (!cedbc) {
610  assert(par_bc != 0x0) ;
611  assert(par_bc->get_n_tbl_mod() >= 2) ;
612  }
613 #endif
614  int nt = mgrid.get_nt(0) ;
615  int np = mgrid.get_np(0) ;
616 
617  assert (tilde_b.get_etat() != ETATNONDEF) ;
618  assert (hh.get_etat() != ETATNONDEF) ;
619 
620  Scalar source = tilde_b ;
621  Scalar source_coq = tilde_b ;
622  source_coq.annule_domain(0) ;
623  if (ced)
624  source_coq.annule_domain(nzm1) ;
625  source_coq.mult_r() ;
626  source.set_spectral_va().ylm() ;
627  source_coq.set_spectral_va().ylm() ;
628  bool bnull = (tilde_b.get_etat() == ETATZERO) ;
629 
630  assert(hh.check_dzpuis(0)) ;
631  Scalar hoverr = hh ;
632  hoverr.div_r_dzpuis(2) ;
633  hoverr.set_spectral_va().ylm() ;
634  Scalar dhdr = hh.dsdr() ;
635  dhdr.set_spectral_va().ylm() ;
636  Scalar h_coq = hh ;
637  h_coq.set_spectral_va().ylm() ;
638  Scalar dh_coq = hh.dsdr() ;
639  dh_coq.mult_r_dzpuis(0) ;
640  dh_coq.set_spectral_va().ylm() ;
641  bool hnull = (hh.get_etat() == ETATZERO) ;
642 
643  Base_val base = (bnull ? hoverr.get_spectral_base() : source.get_spectral_base()) ;
644  base.mult_x() ;
645  int lmax = base.give_lmax(mgrid, 0) + 1;
646 
647  bool need_calculation = true ;
648  if (par_mat != 0x0) {
649  bool param_new = false ;
650  if ((par_mat->get_n_int_mod() >= 4)
651  &&(par_mat->get_n_tbl_mod()>=1)
652  &&(par_mat->get_n_matrice_mod()>=1)
653  &&(par_mat->get_n_itbl_mod()>=1)) {
654  if (par_mat->get_int_mod(0) < nz_bc) param_new = true ;
655  if (par_mat->get_int_mod(1) != lmax) param_new = true ;
656  if (par_mat->get_int_mod(2) != mgrid.get_type_t() ) param_new = true ;
657  if (par_mat->get_int_mod(3) != mgrid.get_type_p() ) param_new = true ;
658  if (par_mat->get_itbl_mod(0)(0) != mgrid.get_nr(0)) param_new = true ;
659  if (fabs(par_mat->get_tbl_mod(0)(0) - mp_aff->get_alpha()[0]) > 2.e-15)
660  param_new = true ;
661  for (int l=1; l<= n_shell; l++) {
662  if (par_mat->get_itbl_mod(0)(l) != mgrid.get_nr(l)) param_new = true ;
663  if (fabs(par_mat->get_tbl_mod(0)(l) - mp_aff->get_beta()[l] /
664  mp_aff->get_alpha()[l]) > 2.e-15) param_new = true ;
665  }
666  if (ced) {
667  if (par_mat->get_itbl_mod(0)(nzm1) != mgrid.get_nr(nzm1)) param_new = true ;
668  if (fabs(par_mat->get_tbl_mod(0)(nzm1) - mp_aff->get_alpha()[nzm1]) > 2.e-15)
669  param_new = true ;
670  }
671  }
672  else{
673  param_new = true ;
674  }
675  if (param_new) {
676  par_mat->clean_all() ;
677  int* nz_bc_new = new int(nz_bc) ;
678  par_mat->add_int_mod(*nz_bc_new, 0) ;
679  int* lmax_new = new int(lmax) ;
680  par_mat->add_int_mod(*lmax_new, 1) ;
681  int* type_t_new = new int(mgrid.get_type_t()) ;
682  par_mat->add_int_mod(*type_t_new, 2) ;
683  int* type_p_new = new int(mgrid.get_type_p()) ;
684  par_mat->add_int_mod(*type_p_new, 3) ;
685  Itbl* pnr = new Itbl(nz) ;
686  pnr->set_etat_qcq() ;
687  par_mat->add_itbl_mod(*pnr) ;
688  for (int l=0; l<nz; l++)
689  pnr->set(l) = mgrid.get_nr(l) ;
690  Tbl* palpha = new Tbl(nz) ;
691  palpha->set_etat_qcq() ;
692  par_mat->add_tbl_mod(*palpha) ;
693  palpha->set(0) = mp_aff->get_alpha()[0] ;
694  for (int l=1; l<nzm1; l++)
695  palpha->set(l) = mp_aff->get_beta()[l] / mp_aff->get_alpha()[l] ;
696  palpha->set(nzm1) = mp_aff->get_alpha()[nzm1] ;
697  }
698  else need_calculation = false ;
699  }
700 
701  hrr.set_etat_qcq() ;
702  hrr.set_spectral_base(base) ;
704  hrr.set_spectral_va().c_cf->annule_hard() ;
705  tilde_eta.annule_hard() ;
706  tilde_eta.set_spectral_base(base) ;
707  tilde_eta.set_spectral_va().set_etat_cf_qcq() ;
708  tilde_eta.set_spectral_va().c_cf->annule_hard() ;
709  ww.annule_hard() ;
710  ww.set_spectral_base(base) ;
712  ww.set_spectral_va().c_cf->annule_hard() ;
713 
714  sol_Dirac_l01(hh, hrr, tilde_eta, par_mat) ;
715  tilde_eta.annule_l(0,0, true) ;
716 
717  Mtbl_cf sol_part_hrr(mgrid, base) ; sol_part_hrr.annule_hard() ;
718  Mtbl_cf sol_part_eta(mgrid, base) ; sol_part_eta.annule_hard() ;
719  Mtbl_cf sol_part_w(mgrid, base) ; sol_part_w.annule_hard() ;
720  Mtbl_cf sol_hom1_hrr(mgrid, base) ; sol_hom1_hrr.annule_hard() ;
721  Mtbl_cf sol_hom1_eta(mgrid, base) ; sol_hom1_eta.annule_hard() ;
722  Mtbl_cf sol_hom1_w(mgrid, base) ; sol_hom1_w.annule_hard() ;
723  Mtbl_cf sol_hom2_hrr(mgrid, base) ; sol_hom2_hrr.annule_hard() ;
724  Mtbl_cf sol_hom2_eta(mgrid, base) ; sol_hom2_eta.annule_hard() ;
725  Mtbl_cf sol_hom2_w(mgrid, base) ; sol_hom2_w.annule_hard() ;
726  Mtbl_cf sol_hom3_hrr(mgrid, base) ; sol_hom3_hrr.annule_hard() ;
727  Mtbl_cf sol_hom3_eta(mgrid, base) ; sol_hom3_eta.annule_hard() ;
728  Mtbl_cf sol_hom3_w(mgrid, base) ; sol_hom3_w.annule_hard() ;
729 
730  int l_q, m_q, base_r ;
731  Itbl mat_done(lmax) ;
732 
733  //---------------
734  //-- NUCLEUS ---
735  //---------------
736  {int lz = 0 ;
737  int nr = mgrid.get_nr(lz) ;
738  double alpha = mp_aff->get_alpha()[lz] ;
739  Matrice ope(3*nr, 3*nr) ;
740  int ind2 = 2*nr ;
741  if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
742 
743  for (int k=0 ; k<np+1 ; k++) {
744  for (int j=0 ; j<nt ; j++) {
745  // quantic numbers and spectral bases
746  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
747  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
748  if (need_calculation) {
749  ope.set_etat_qcq() ;
750  Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
751  Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
752 
753  for (int lin=0; lin<nr; lin++)
754  for (int col=0; col<nr; col++)
755  ope.set(lin,col) = md(lin,col) + 3*ms(lin,col) ;
756  for (int lin=0; lin<nr; lin++)
757  for (int col=0; col<nr; col++)
758  ope.set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
759  for (int lin=0; lin<nr; lin++)
760  for (int col=0; col<nr; col++)
761  ope.set(lin,col+2*nr) = 0 ;
762  for (int lin=0; lin<nr; lin++)
763  for (int col=0; col<nr; col++)
764  ope.set(lin+nr,col) = -0.5*ms(lin,col) ;
765  for (int lin=0; lin<nr; lin++)
766  for (int col=0; col<nr; col++)
767  ope.set(lin+nr,col+nr) = md(lin,col) + 3*ms(lin,col) ;
768  for (int lin=0; lin<nr; lin++)
769  for (int col=0; col<nr; col++)
770  ope.set(lin+nr,col+2*nr) = (2. - l_q*(l_q+1))*ms(lin,col) ;
771  for (int lin=0; lin<nr; lin++)
772  for (int col=0; col<nr; col++)
773  ope.set(lin+2*nr,col) = -0.5*md(lin,col)/double(l_q+1)
774  - 0.5*double(l_q+4)/double(l_q+1)*ms(lin,col) ;
775  for (int lin=0; lin<nr; lin++)
776  for (int col=0; col<nr; col++)
777  ope.set(lin+2*nr,col+nr) = -2*ms(lin,col) ;
778  for (int lin=0; lin<nr; lin++)
779  for (int col=0; col<nr; col++)
780  ope.set(lin+2*nr,col+2*nr) = (l_q+2)*md(lin,col)
781  + l_q*(l_q+2)*ms(lin,col) ;
782 
783  ope *= 1./alpha ;
784  for (int col=0; col<3*nr; col++)
785  if (l_q>2) ope.set(ind2+nr-2, col) = 0 ;
786  for (int col=0; col<3*nr; col++) {
787  ope.set(nr-1, col) = 0 ;
788  ope.set(2*nr-1, col) = 0 ;
789  ope.set(3*nr-1, col) = 0 ;
790  }
791  int pari = 1 ;
792  if (base_r == R_CHEBP) {
793  for (int col=0; col<nr; col++) {
794  ope.set(nr-1, col) = pari ;
795  ope.set(2*nr-1, col+nr) = pari ;
796  ope.set(3*nr-1, col+2*nr) = pari ;
797  pari = - pari ;
798  }
799  }
800  else { //In the odd case, the last coefficient must be zero!
801  ope.set(nr-1, nr-1) = 1 ;
802  ope.set(2*nr-1, 2*nr-1) = 1 ;
803  ope.set(3*nr-1, 3*nr-1) = 1 ;
804  }
805  if (l_q>2)
806  ope.set(ind2+nr-2, ind2) = 1 ;
807 
808  ope.set_lu() ;
809  if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
810  Matrice* pope = new Matrice(ope) ;
811  par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
812  mat_done.set(l_q) = 1 ;
813  }
814  } //End of case when a calculation is needed
815 
816  const Matrice& oper = (par_mat == 0x0 ? ope :
817  par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
818  Tbl sec(3*nr) ;
819  sec.set_etat_qcq() ;
820  if (hnull) {
821  for (int lin=0; lin<2*nr; lin++)
822  sec.set(lin) = 0 ;
823  for (int lin=0; lin<nr; lin++)
824  sec.set(2*nr+lin) = (*source.get_spectral_va().c_cf)
825  (lz, k, j, lin) ;
826  }
827  else {
828  for (int lin=0; lin<nr; lin++)
829  sec.set(lin) = (*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) ;
830  for (int lin=0; lin<nr; lin++)
831  sec.set(lin+nr) = -0.5*(*hoverr.get_spectral_va().c_cf)
832  (lz, k, j, lin) ;
833  if (bnull) {
834  for (int lin=0; lin<nr; lin++)
835  sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
836  (*dhdr.get_spectral_va().c_cf)(lz, k, j, lin)
837  + (l_q+2)*(*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) ) ;
838  }
839  else {
840  for (int lin=0; lin<nr; lin++)
841  sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
842  (*dhdr.get_spectral_va().c_cf)(lz, k, j, lin)
843  + (l_q+2)*(*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) )
844  + (*source.get_spectral_va().c_cf)(lz, k, j, lin) ;
845  }
846  }
847  if (l_q>2) sec.set(ind2+nr-2) = 0 ;
848  sec.set(3*nr-1) = 0 ;
849  Tbl sol = oper.inverse(sec) ;
850  for (int i=0; i<nr; i++) {
851  sol_part_hrr.set(lz, k, j, i) = sol(i) ;
852  sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
853  sol_part_w.set(lz, k, j, i) = sol(i+2*nr) ;
854  }
855  sec.annule_hard() ;
856  if (l_q>2) {
857  sec.set(ind2+nr-2) = 1 ;
858  sol = oper.inverse(sec) ;
859  }
860  else { //Homogeneous solution put in by hand in the case l=2
861  sol.annule_hard() ;
862  sol.set(0) = 4 ;
863  sol.set(nr) = 2 ;
864  sol.set(2*nr) = 1 ;
865  }
866  for (int i=0; i<nr; i++) {
867  sol_hom3_hrr.set(lz, k, j, i) = sol(i) ;
868  sol_hom3_eta.set(lz, k, j, i) = sol(i+nr) ;
869  sol_hom3_w.set(lz, k, j, i) = sol(i+2*nr) ;
870  }
871  }
872  }
873  }
874  }
875 
876  //-------------
877  // -- Shells --
878  //-------------
879 
880  for (int lz=1; lz<= n_shell; lz++) {
881  if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
882  int nr = mgrid.get_nr(lz) ;
883  int ind0 = 0 ;
884  int ind1 = nr ;
885  int ind2 = 2*nr ;
886  double alpha = mp_aff->get_alpha()[lz] ;
887  double ech = mp_aff->get_beta()[lz] / alpha ;
888  Matrice ope(3*nr, 3*nr) ;
889 
890  for (int k=0 ; k<np+1 ; k++) {
891  for (int j=0 ; j<nt ; j++) {
892  // quantic numbers and spectral bases
893  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
894  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
895  if (need_calculation) {
896  ope.set_etat_qcq() ;
897  Diff_xdsdx oxd(base_r, nr) ; const Matrice& mxd = oxd.get_matrice() ;
898  Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
899  Diff_id oid(base_r, nr) ; const Matrice& mid = oid.get_matrice() ;
900 
901  for (int lin=0; lin<nr; lin++)
902  for (int col=0; col<nr; col++)
903  ope.set(lin,col) = mxd(lin,col) + ech*md(lin,col)
904  + 3*mid(lin,col) ;
905  for (int lin=0; lin<nr; lin++)
906  for (int col=0; col<nr; col++)
907  ope.set(lin,col+nr) = -l_q*(l_q+1)*mid(lin,col) ;
908  for (int lin=0; lin<nr; lin++)
909  for (int col=0; col<nr; col++)
910  ope.set(lin,col+2*nr) = 0 ;
911  for (int lin=0; lin<nr; lin++)
912  for (int col=0; col<nr; col++)
913  ope.set(lin+nr,col) = -0.5*mid(lin,col) ;
914  for (int lin=0; lin<nr; lin++)
915  for (int col=0; col<nr; col++)
916  ope.set(lin+nr,col+nr) = mxd(lin,col) + ech*md(lin,col)
917  + 3*mid(lin,col) ;
918  for (int lin=0; lin<nr; lin++)
919  for (int col=0; col<nr; col++)
920  ope.set(lin+nr,col+2*nr) = (2. - l_q*(l_q+1))*mid(lin,col) ;
921  for (int lin=0; lin<nr; lin++)
922  for (int col=0; col<nr; col++)
923  ope.set(lin+2*nr,col) =
924  -0.5/double(l_q+1)*(mxd(lin,col) + ech*md(lin,col)
925  + double(l_q+4)*mid(lin,col)) ;
926  for (int lin=0; lin<nr; lin++)
927  for (int col=0; col<nr; col++)
928  ope.set(lin+2*nr,col+nr) = -2*mid(lin,col) ;
929  for (int lin=0; lin<nr; lin++)
930  for (int col=0; col<nr; col++)
931  ope.set(lin+2*nr,col+2*nr) =
932  double(l_q+2)*(mxd(lin,col) + ech*md(lin,col)
933  + l_q*mid(lin,col)) ;
934  for (int col=0; col<3*nr; col++) {
935  ope.set(ind0+nr-1, col) = 0 ;
936  ope.set(ind1+nr-1, col) = 0 ;
937  ope.set(ind2+nr-1, col) = 0 ;
938  }
939  ope.set(ind0+nr-1, ind0) = 1 ;
940  ope.set(ind1+nr-1, ind1) = 1 ;
941  ope.set(ind2+nr-1, ind2) = 1 ;
942 
943  ope.set_lu() ;
944  if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
945  Matrice* pope = new Matrice(ope) ;
946  par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
947  mat_done.set(l_q) = 1 ;
948  }
949  } //End of case when a calculation is needed
950  const Matrice& oper = (par_mat == 0x0 ? ope :
951  par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
952  Tbl sec(3*nr) ;
953  sec.set_etat_qcq() ;
954  if (hnull) {
955  for (int lin=0; lin<2*nr; lin++)
956  sec.set(lin) = 0 ;
957  for (int lin=0; lin<nr; lin++)
958  sec.set(2*nr+lin) = (*source_coq.get_spectral_va().c_cf)
959  (lz, k, j, lin) ;
960  }
961  else {
962  for (int lin=0; lin<nr; lin++)
963  sec.set(lin) = (*h_coq.get_spectral_va().c_cf)(lz, k, j, lin) ;
964  for (int lin=0; lin<nr; lin++)
965  sec.set(lin+nr) = -0.5*(*h_coq.get_spectral_va().c_cf)
966  (lz, k, j, lin) ;
967  if (bnull) {
968  for (int lin=0; lin<nr; lin++)
969  sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
970  (*dh_coq.get_spectral_va().c_cf)(lz, k, j, lin)
971  + (l_q+2)*(*h_coq.get_spectral_va().c_cf)(lz, k, j, lin) ) ;
972  }
973  else {
974  for (int lin=0; lin<nr; lin++)
975  sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
976  (*dh_coq.get_spectral_va().c_cf)(lz, k, j, lin)
977  + (l_q+2)*(*h_coq.get_spectral_va().c_cf)(lz, k, j, lin) )
978  + (*source_coq.get_spectral_va().c_cf)(lz, k, j, lin) ;
979  }
980  }
981  sec.set(ind0+nr-1) = 0 ;
982  sec.set(ind1+nr-1) = 0 ;
983  sec.set(ind2+nr-1) = 0 ;
984  Tbl sol = oper.inverse(sec) ;
985  for (int i=0; i<nr; i++) {
986  sol_part_hrr.set(lz, k, j, i) = sol(i) ;
987  sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
988  sol_part_w.set(lz, k, j, i) = sol(i+2*nr) ;
989  }
990  sec.annule_hard() ;
991  sec.set(ind0+nr-1) = 1 ;
992  sol = oper.inverse(sec) ;
993  for (int i=0; i<nr; i++) {
994  sol_hom1_hrr.set(lz, k, j, i) = sol(i) ;
995  sol_hom1_eta.set(lz, k, j, i) = sol(i+nr) ;
996  sol_hom1_w.set(lz, k, j, i) = sol(i+2*nr) ;
997  }
998  sec.set(ind0+nr-1) = 0 ;
999  sec.set(ind1+nr-1) = 1 ;
1000  sol = oper.inverse(sec) ;
1001  for (int i=0; i<nr; i++) {
1002  sol_hom2_hrr.set(lz, k, j, i) = sol(i) ;
1003  sol_hom2_eta.set(lz, k, j, i) = sol(i+nr) ;
1004  sol_hom2_w.set(lz, k, j, i) = sol(i+2*nr) ;
1005  }
1006  sec.set(ind1+nr-1) = 0 ;
1007  sec.set(ind2+nr-1) = 1 ;
1008  sol = oper.inverse(sec) ;
1009  for (int i=0; i<nr; i++) {
1010  sol_hom3_hrr.set(lz, k, j, i) = sol(i) ;
1011  sol_hom3_eta.set(lz, k, j, i) = sol(i+nr) ;
1012  sol_hom3_w.set(lz, k, j, i) = sol(i+2*nr) ;
1013  }
1014  }
1015  }
1016  }
1017  }
1018 
1019  //------------------------------
1020  // Compactified external domain
1021  //------------------------------
1022  if (cedbc) {int lz = nzm1 ;
1023  if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
1024  int nr = mgrid.get_nr(lz) ;
1025  int ind0 = 0 ;
1026  int ind1 = nr ;
1027  int ind2 = 2*nr ;
1028  double alpha = mp_aff->get_alpha()[lz] ;
1029  Matrice ope(3*nr, 3*nr) ;
1030 
1031  for (int k=0 ; k<np+1 ; k++) {
1032  for (int j=0 ; j<nt ; j++) {
1033  // quantic numbers and spectral bases
1034  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
1035  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
1036  if (need_calculation) {
1037  ope.set_etat_qcq() ;
1038  Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
1039  Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
1040 
1041  for (int lin=0; lin<nr; lin++)
1042  for (int col=0; col<nr; col++)
1043  ope.set(lin,col) = - md(lin,col) + 3*ms(lin,col) ;
1044  for (int lin=0; lin<nr; lin++)
1045  for (int col=0; col<nr; col++)
1046  ope.set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
1047  for (int lin=0; lin<nr; lin++)
1048  for (int col=0; col<nr; col++)
1049  ope.set(lin,col+2*nr) = 0 ;
1050  for (int lin=0; lin<nr; lin++)
1051  for (int col=0; col<nr; col++)
1052  ope.set(lin+nr,col) = -0.5*ms(lin,col) ;
1053  for (int lin=0; lin<nr; lin++)
1054  for (int col=0; col<nr; col++)
1055  ope.set(lin+nr,col+nr) = -md(lin,col) + 3*ms(lin,col) ;
1056  for (int lin=0; lin<nr; lin++)
1057  for (int col=0; col<nr; col++)
1058  ope.set(lin+nr,col+2*nr) = (2. - l_q*(l_q+1))*ms(lin,col) ;
1059  for (int lin=0; lin<nr; lin++)
1060  for (int col=0; col<nr; col++)
1061  ope.set(lin+2*nr,col) = 0.5*md(lin,col)/double(l_q+1)
1062  - 0.5*double(l_q+4)/double(l_q+1)*ms(lin,col) ;
1063  for (int lin=0; lin<nr; lin++)
1064  for (int col=0; col<nr; col++)
1065  ope.set(lin+2*nr,col+nr) = -2*ms(lin,col) ;
1066  for (int lin=0; lin<nr; lin++)
1067  for (int col=0; col<nr; col++)
1068  ope.set(lin+2*nr,col+2*nr) = -(l_q+2)*md(lin,col)
1069  + l_q*(l_q+2)*ms(lin,col) ;
1070  ope *= 1./alpha ;
1071  for (int col=0; col<3*nr; col++) {
1072  ope.set(ind0+nr-2, col) = 0 ;
1073  ope.set(ind0+nr-1, col) = 0 ;
1074  ope.set(ind1+nr-2, col) = 0 ;
1075  ope.set(ind1+nr-1, col) = 0 ;
1076  ope.set(ind2+nr-1, col) = 0 ;
1077  }
1078  for (int col=0; col<nr; col++) {
1079  ope.set(ind0+nr-1, col+ind0) = 1 ;
1080  ope.set(ind1+nr-1, col+ind1) = 1 ;
1081  ope.set(ind2+nr-1, col+ind2) = 1 ;
1082  }
1083  ope.set(ind0+nr-2, ind0+1) = 1 ;
1084  ope.set(ind1+nr-2, ind1+2) = 1 ;
1085 
1086  ope.set_lu() ;
1087  if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
1088  Matrice* pope = new Matrice(ope) ;
1089  par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
1090  mat_done.set(l_q) = 1 ;
1091  }
1092  } //End of case when a calculation is needed
1093  const Matrice& oper = (par_mat == 0x0 ? ope :
1094  par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
1095 
1096  Tbl sec(3*nr) ;
1097  sec.set_etat_qcq() ;
1098  if (hnull) {
1099  for (int lin=0; lin<2*nr; lin++)
1100  sec.set(lin) = 0 ;
1101  for (int lin=0; lin<nr; lin++)
1102  sec.set(2*nr+lin) = (*source.get_spectral_va().c_cf)
1103  (lz, k, j, lin) ;
1104  }
1105  else {
1106  for (int lin=0; lin<nr; lin++)
1107  sec.set(lin) = (*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) ;
1108  for (int lin=0; lin<nr; lin++)
1109  sec.set(lin+nr) = -0.5*(*hoverr.get_spectral_va().c_cf)
1110  (lz, k, j, lin) ;
1111  if (bnull) {
1112  for (int lin=0; lin<nr; lin++)
1113  sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
1114  (*dhdr.get_spectral_va().c_cf)(lz, k, j, lin)
1115  + (l_q+2)*(*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) ) ;
1116  }
1117  else {
1118  for (int lin=0; lin<nr; lin++)
1119  sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
1120  (*dhdr.get_spectral_va().c_cf)(lz, k, j, lin)
1121  + (l_q+2)*(*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) )
1122  + (*source.get_spectral_va().c_cf)(lz, k, j, lin) ;
1123  }
1124  }
1125  sec.set(ind0+nr-2) = 0 ;
1126  sec.set(ind0+nr-1) = 0 ;
1127  sec.set(ind1+nr-1) = 0 ;
1128  sec.set(ind1+nr-2) = 0 ;
1129  sec.set(ind2+nr-1) = 0 ;
1130  Tbl sol = oper.inverse(sec) ;
1131  for (int i=0; i<nr; i++) {
1132  sol_part_hrr.set(lz, k, j, i) = sol(i) ;
1133  sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
1134  sol_part_w.set(lz, k, j, i) = sol(i+2*nr) ;
1135  }
1136  sec.annule_hard() ;
1137  sec.set(ind0+nr-2) = 1 ;
1138  sol = oper.inverse(sec) ;
1139  for (int i=0; i<nr; i++) {
1140  sol_hom1_hrr.set(lz, k, j, i) = sol(i) ;
1141  sol_hom1_eta.set(lz, k, j, i) = sol(i+nr) ;
1142  sol_hom1_w.set(lz, k, j, i) = sol(i+2*nr) ;
1143  }
1144  sec.set(ind0+nr-2) = 0 ;
1145  sec.set(ind1+nr-2) = 1 ;
1146  sol = oper.inverse(sec) ;
1147  for (int i=0; i<nr; i++) {
1148  sol_hom2_hrr.set(lz, k, j, i) = sol(i) ;
1149  sol_hom2_eta.set(lz, k, j, i) = sol(i+nr) ;
1150  sol_hom2_w.set(lz, k, j, i) = sol(i+2*nr) ;
1151  }
1152  }
1153  }
1154  }
1155  }
1156 
1157  int taille = 3*nz_bc + 1 ;
1158  if (cedbc) taille-- ;
1159  Mtbl_cf& mhrr = *hrr.set_spectral_va().c_cf ;
1160  Mtbl_cf& meta = *tilde_eta.set_spectral_va().c_cf ;
1161  Mtbl_cf& mw = *ww.set_spectral_va().c_cf ;
1162 
1163  Tbl sec_membre(taille) ;
1164  Matrice systeme(taille, taille) ;
1165  int ligne ; int colonne ;
1166  Tbl pipo(1) ;
1167  const Tbl& hrrb = (cedbc ? pipo : par_bc->get_tbl_mod(1) );
1168  double chrr = (cedbc ? 0 : par_bc->get_tbl_mod()(4) ) ;
1169  double dhrr = (cedbc ? 0 : par_bc->get_tbl_mod()(5) ) ;
1170  double ceta = (cedbc ? 0 : par_bc->get_tbl_mod()(6) ) ;
1171  double deta = (cedbc ? 0 : par_bc->get_tbl_mod()(7) ) ;
1172  double cw = (cedbc ? 0 : par_bc->get_tbl_mod()(8) ) ;
1173  double dw = (cedbc ? 0 : par_bc->get_tbl_mod()(9) ) ;
1174  Mtbl_cf dhom1_hrr = sol_hom1_hrr ;
1175  Mtbl_cf dhom2_hrr = sol_hom2_hrr ;
1176  Mtbl_cf dhom3_hrr = sol_hom3_hrr ;
1177  Mtbl_cf dpart_hrr = sol_part_hrr ;
1178  Mtbl_cf dhom1_eta = sol_hom1_eta ;
1179  Mtbl_cf dhom2_eta = sol_hom2_eta ;
1180  Mtbl_cf dhom3_eta = sol_hom3_eta ;
1181  Mtbl_cf dpart_eta = sol_part_eta ;
1182  Mtbl_cf dhom1_w = sol_hom1_w ;
1183  Mtbl_cf dhom2_w = sol_hom2_w ;
1184  Mtbl_cf dhom3_w = sol_hom3_w ;
1185  Mtbl_cf dpart_w = sol_part_w ;
1186  if (!cedbc) {
1187  dhom1_hrr.dsdx() ; dhom1_eta.dsdx() ; dhom1_w.dsdx() ;
1188  dhom2_hrr.dsdx() ; dhom2_eta.dsdx() ; dhom2_w.dsdx() ;
1189  dhom3_hrr.dsdx() ; dhom3_eta.dsdx() ; dhom3_w.dsdx() ;
1190  dpart_hrr.dsdx() ; dpart_eta.dsdx() ; dpart_w.dsdx() ;
1191  }
1192  // Loop on l and m
1193  //----------------
1194  for (int k=0 ; k<np+1 ; k++)
1195  for (int j=0 ; j<nt ; j++) {
1196  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
1197  if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
1198  ligne = 0 ;
1199  colonne = 0 ;
1200  systeme.annule_hard() ;
1201  sec_membre.annule_hard() ;
1202 
1203  //Nucleus
1204  int nr = mgrid.get_nr(0) ;
1205  if (nz_bc >0 ) {
1206  systeme.set(ligne, colonne) = sol_hom3_hrr.val_out_bound_jk(0, j, k) ;
1207  sec_membre.set(ligne) = -sol_part_hrr.val_out_bound_jk(0, j, k) ;
1208  ligne++ ;
1209 
1210  systeme.set(ligne, colonne) = sol_hom3_eta.val_out_bound_jk(0, j, k) ;
1211  sec_membre.set(ligne) = -sol_part_eta.val_out_bound_jk(0, j, k) ;
1212  ligne++ ;
1213 
1214  systeme.set(ligne, colonne) = sol_hom3_w.val_out_bound_jk(0, j, k) ;
1215  sec_membre.set(ligne) = -sol_part_w.val_out_bound_jk(0, j, k) ;
1216  colonne++ ;
1217  }
1218  //shells
1219  for (int zone=1 ; zone<nz_bc ; zone++) {
1220  nr = mgrid.get_nr(zone) ;
1221  ligne -= 2 ;
1222 
1223  //Condition at x = -1
1224  systeme.set(ligne, colonne) =
1225  - sol_hom1_hrr.val_in_bound_jk(zone, j, k) ;
1226  systeme.set(ligne, colonne+1) =
1227  - sol_hom2_hrr.val_in_bound_jk(zone, j, k) ;
1228  systeme.set(ligne, colonne+2) =
1229  - sol_hom3_hrr.val_in_bound_jk(zone, j, k) ;
1230 
1231  sec_membre.set(ligne) += sol_part_hrr.val_in_bound_jk(zone, j, k) ;
1232  ligne++ ;
1233 
1234  systeme.set(ligne, colonne) =
1235  - sol_hom1_eta.val_in_bound_jk(zone, j, k) ;
1236  systeme.set(ligne, colonne+1) =
1237  - sol_hom2_eta.val_in_bound_jk(zone, j, k) ;
1238  systeme.set(ligne, colonne+2) =
1239  - sol_hom3_eta.val_in_bound_jk(zone, j, k) ;
1240 
1241  sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(zone, j, k) ;
1242  ligne++ ;
1243 
1244  systeme.set(ligne, colonne) =
1245  - sol_hom1_w.val_in_bound_jk(zone, j, k) ;
1246  systeme.set(ligne, colonne+1) =
1247  - sol_hom2_w.val_in_bound_jk(zone, j, k) ;
1248  systeme.set(ligne, colonne+2) =
1249  - sol_hom3_w.val_in_bound_jk(zone, j, k) ;
1250 
1251  sec_membre.set(ligne) += sol_part_w.val_in_bound_jk(zone, j, k) ;
1252  ligne++ ;
1253 
1254  // Condition at x=1
1255  systeme.set(ligne, colonne) =
1256  sol_hom1_hrr.val_out_bound_jk(zone, j, k) ;
1257  systeme.set(ligne, colonne+1) =
1258  sol_hom2_hrr.val_out_bound_jk(zone, j, k) ;
1259  systeme.set(ligne, colonne+2) =
1260  sol_hom3_hrr.val_out_bound_jk(zone, j, k) ;
1261 
1262  sec_membre.set(ligne) -= sol_part_hrr.val_out_bound_jk(zone, j, k) ;
1263  ligne++ ;
1264 
1265  systeme.set(ligne, colonne) =
1266  sol_hom1_eta.val_out_bound_jk(zone, j, k) ;
1267  systeme.set(ligne, colonne+1) =
1268  sol_hom2_eta.val_out_bound_jk(zone, j, k) ;
1269  systeme.set(ligne, colonne+2) =
1270  sol_hom3_eta.val_out_bound_jk(zone, j, k) ;
1271 
1272  sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone, j, k) ;
1273  ligne++ ;
1274 
1275  systeme.set(ligne, colonne) =
1276  sol_hom1_w.val_out_bound_jk(zone, j, k) ;
1277  systeme.set(ligne, colonne+1) =
1278  sol_hom2_w.val_out_bound_jk(zone, j, k) ;
1279  systeme.set(ligne, colonne+2) =
1280  sol_hom3_w.val_out_bound_jk(zone, j, k) ;
1281 
1282  sec_membre.set(ligne) -= sol_part_w.val_out_bound_jk(zone, j, k) ;
1283 
1284  colonne += 3 ;
1285  }
1286 
1287  //Last domain
1288  nr = mgrid.get_nr(nz_bc) ;
1289  double alpha = mp_aff->get_alpha()[nz_bc] ;
1290  if (nz_bc>0) {
1291  ligne -= 2 ;
1292 
1293  //Condition at x = -1
1294  systeme.set(ligne, colonne) =
1295  - sol_hom1_hrr.val_in_bound_jk(nz_bc, j, k) ;
1296  systeme.set(ligne, colonne+1) =
1297  - sol_hom2_hrr.val_in_bound_jk(nz_bc, j, k) ;
1298  if (!cedbc) systeme.set(ligne, colonne+2) =
1299  - sol_hom3_hrr.val_in_bound_jk(nz_bc, j, k) ;
1300 
1301  sec_membre.set(ligne) += sol_part_hrr.val_in_bound_jk(nz_bc, j, k) ;
1302  ligne++ ;
1303 
1304  systeme.set(ligne, colonne) =
1305  - sol_hom1_eta.val_in_bound_jk(nz_bc, j, k) ;
1306  systeme.set(ligne, colonne+1) =
1307  - sol_hom2_eta.val_in_bound_jk(nz_bc, j, k) ;
1308  if (!cedbc) systeme.set(ligne, colonne+2) =
1309  - sol_hom3_eta.val_in_bound_jk(nz_bc, j, k) ;
1310 
1311  sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(nz_bc, j, k) ;
1312  ligne++ ;
1313 
1314  systeme.set(ligne, colonne) =
1315  - sol_hom1_w.val_in_bound_jk(nz_bc, j, k) ;
1316  systeme.set(ligne, colonne+1) =
1317  - sol_hom2_w.val_in_bound_jk(nz_bc, j, k) ;
1318  if (!cedbc) systeme.set(ligne, colonne+2) =
1319  - sol_hom3_w.val_in_bound_jk(nz_bc, j, k) ;
1320 
1321  sec_membre.set(ligne) += sol_part_w.val_in_bound_jk(nz_bc, j, k) ;
1322  ligne++ ;
1323  }
1324  if (!cedbc) {//Special condition at x=1
1325  if (nz_bc > 0) {
1326  systeme.set(ligne, colonne) =
1327  chrr*sol_hom1_hrr.val_out_bound_jk(nz_bc, j, k)
1328  + dhrr*dhom1_hrr.val_out_bound_jk(nz_bc, j, k) / alpha
1329  + ceta*sol_hom1_eta.val_out_bound_jk(nz_bc, j, k)
1330  + deta*dhom1_eta.val_out_bound_jk(nz_bc, j, k) / alpha
1331  + cw*sol_hom1_w.val_out_bound_jk(nz_bc, j, k)
1332  + dw*dhom1_w.val_out_bound_jk(nz_bc, j, k) / alpha ;
1333  systeme.set(ligne, colonne+1) =
1334  chrr*sol_hom2_hrr.val_out_bound_jk(nz_bc, j, k)
1335  + dhrr*dhom2_hrr.val_out_bound_jk(nz_bc, j, k) / alpha
1336  + ceta*sol_hom2_eta.val_out_bound_jk(nz_bc, j, k)
1337  + deta*dhom2_eta.val_out_bound_jk(nz_bc, j, k) / alpha
1338  + cw*sol_hom2_w.val_out_bound_jk(nz_bc, j, k)
1339  + dw*dhom2_w.val_out_bound_jk(nz_bc, j, k) / alpha ;
1340  }
1341  else { // Only nucleus
1342  assert(ligne == 0) ;
1343  colonne = -2 ;
1344  }
1345  systeme.set(ligne, colonne+2) =
1346  chrr*sol_hom3_hrr.val_out_bound_jk(nz_bc, j, k)
1347  + dhrr*dhom3_hrr.val_out_bound_jk(nz_bc, j, k) / alpha
1348  + ceta*sol_hom3_eta.val_out_bound_jk(nz_bc, j, k)
1349  + deta*dhom3_eta.val_out_bound_jk(nz_bc, j, k) / alpha
1350  + cw*sol_hom3_w.val_out_bound_jk(nz_bc, j, k)
1351  + dw*dhom3_w.val_out_bound_jk(nz_bc, j, k) / alpha ;
1352 
1353  sec_membre.set(ligne) -= chrr*sol_part_hrr.val_out_bound_jk(nz_bc, j, k)
1354  + dhrr*dpart_hrr.val_out_bound_jk(nz_bc, j, k)/alpha
1355  + ceta*sol_part_eta.val_out_bound_jk(nz_bc, j, k)
1356  + deta*dpart_eta.val_out_bound_jk(nz_bc, j, k)/alpha
1357  + cw*sol_part_w.val_out_bound_jk(nz_bc, j, k)
1358  + dw*dpart_w.val_out_bound_jk(nz_bc, j, k)/alpha
1359  - hrrb(k, j) ;
1360  }
1361 
1362  // Solution of the system giving the coefficients for the homogeneous
1363  // solutions
1364  //-------------------------------------------------------------------
1365  systeme.set_lu() ;
1366  Tbl facteur = systeme.inverse(sec_membre) ;
1367  int conte = 0 ;
1368 
1369  // everything is put to the right place,
1370  //---------------------------------------
1371  nr = mgrid.get_nr(0) ; //nucleus
1372  for (int i=0 ; i<nr ; i++) {
1373  mhrr.set(0, k, j, i) = sol_part_hrr(0, k, j, i)
1374  + facteur(conte)*sol_hom3_hrr(0, k, j, i) ;
1375  meta.set(0, k, j, i) = sol_part_eta(0, k, j, i)
1376  + facteur(conte)*sol_hom3_eta(0, k, j, i) ;
1377  mw.set(0, k, j, i) = sol_part_w(0, k, j, i)
1378  + facteur(conte)*sol_hom3_w(0, k, j, i) ;
1379  }
1380  conte++ ;
1381  for (int zone=1 ; zone<=n_shell ; zone++) { //shells
1382  nr = mgrid.get_nr(zone) ;
1383  for (int i=0 ; i<nr ; i++) {
1384  mhrr.set(zone, k, j, i) = sol_part_hrr(zone, k, j, i)
1385  + facteur(conte)*sol_hom1_hrr(zone, k, j, i)
1386  + facteur(conte+1)*sol_hom2_hrr(zone, k, j, i)
1387  + facteur(conte+2)*sol_hom3_hrr(zone, k, j, i) ;
1388 
1389  meta.set(zone, k, j, i) = sol_part_eta(zone, k, j, i)
1390  + facteur(conte)*sol_hom1_eta(zone, k, j, i)
1391  + facteur(conte+1)*sol_hom2_eta(zone, k, j, i)
1392  + facteur(conte+2)*sol_hom3_eta(zone, k, j, i) ;
1393 
1394  mw.set(zone, k, j, i) = sol_part_w(zone, k, j, i)
1395  + facteur(conte)*sol_hom1_w(zone, k, j, i)
1396  + facteur(conte+1)*sol_hom2_w(zone, k, j, i)
1397  + facteur(conte+2)*sol_hom3_w(zone, k, j, i) ;
1398  }
1399  conte+=3 ;
1400  }
1401  if (cedbc) {
1402  nr = mgrid.get_nr(nzm1) ; //compactified external domain
1403  for (int i=0 ; i<nr ; i++) {
1404  mhrr.set(nzm1, k, j, i) = sol_part_hrr(nzm1, k, j, i)
1405  + facteur(conte)*sol_hom1_hrr(nzm1, k, j, i)
1406  + facteur(conte+1)*sol_hom2_hrr(nzm1, k, j, i) ;
1407 
1408  meta.set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
1409  + facteur(conte)*sol_hom1_eta(nzm1, k, j, i)
1410  + facteur(conte+1)*sol_hom2_eta(nzm1, k, j, i) ;
1411 
1412  mw.set(nzm1, k, j, i) = sol_part_w(nzm1, k, j, i)
1413  + facteur(conte)*sol_hom1_w(nzm1, k, j, i)
1414  + facteur(conte+1)*sol_hom2_w(nzm1, k, j, i) ;
1415  }
1416  }
1417  } // End of nullite_plm
1418  } //End of loop on theta
1419 
1420 
1421  if (hrr.set_spectral_va().c != 0x0)
1422  delete hrr.set_spectral_va().c ;
1423  hrr.set_spectral_va().c = 0x0 ;
1424  hrr.set_spectral_va().ylm_i() ;
1425 
1426  if (tilde_eta.set_spectral_va().c != 0x0)
1427  delete tilde_eta.set_spectral_va().c ;
1428  tilde_eta.set_spectral_va().c = 0x0 ;
1429  tilde_eta.set_spectral_va().ylm_i() ;
1430 
1431  if (ww.set_spectral_va().c != 0x0)
1432  delete ww.set_spectral_va().c ;
1433  ww.set_spectral_va().c = 0x0 ;
1434  ww.set_spectral_va().ylm_i() ;
1435 
1436 }
1437 
1438 void Sym_tensor_trans::sol_Dirac_l01(const Scalar& hh, Scalar& hrr, Scalar& tilde_eta,
1439  Param* par_mat) const {
1440 
1441  const Map_af* mp_aff = dynamic_cast<const Map_af*>(mp) ;
1442  assert(mp_aff != 0x0) ; //Only affine mapping for the moment
1443 
1444  if (hh.get_etat() == ETATZERO) {
1445  hrr.annule_hard() ;
1446  tilde_eta.annule_hard() ;
1447  return ;
1448  }
1449 
1450  const Mg3d& mgrid = *mp_aff->get_mg() ;
1451  int nz = mgrid.get_nzone() ;
1452  assert(mgrid.get_type_r(0) == RARE) ;
1453  assert(mgrid.get_type_r(nz-1) == UNSURR) ;
1454 
1455  int nt = mgrid.get_nt(0) ;
1456  int np = mgrid.get_np(0) ;
1457 
1458  Scalar source = hh ;
1459  source.div_r_dzpuis(2) ;
1460  Scalar source_coq = hh ;
1461  source.set_spectral_va().ylm() ;
1462  source_coq.set_spectral_va().ylm() ;
1463  Base_val base = source.get_spectral_base() ;
1464  base.mult_x() ;
1465  int lmax = base.give_lmax(mgrid, 0) + 1;
1466 
1467  assert (hrr.get_spectral_base() == base) ;
1468  assert (tilde_eta.get_spectral_base() == base) ;
1469  assert (hrr.get_spectral_va().c_cf != 0x0) ;
1470  assert (tilde_eta.get_spectral_va().c_cf != 0x0) ;
1471 
1472  Mtbl_cf sol_part_hrr(mgrid, base) ; sol_part_hrr.annule_hard() ;
1473  Mtbl_cf sol_part_eta(mgrid, base) ; sol_part_eta.annule_hard() ;
1474  Mtbl_cf sol_hom1_hrr(mgrid, base) ; sol_hom1_hrr.annule_hard() ;
1475  Mtbl_cf sol_hom1_eta(mgrid, base) ; sol_hom1_eta.annule_hard() ;
1476  Mtbl_cf sol_hom2_hrr(mgrid, base) ; sol_hom2_hrr.annule_hard() ;
1477  Mtbl_cf sol_hom2_eta(mgrid, base) ; sol_hom2_eta.annule_hard() ;
1478 
1479  bool need_calculation = true ;
1480  if (par_mat != 0x0)
1481  if (par_mat->get_n_matrice_mod() > 0)
1482  if (&par_mat->get_matrice_mod(0) != 0x0) need_calculation = false ;
1483 
1484  int l_q, m_q, base_r ;
1485  Itbl mat_done(lmax) ;
1486 
1487  //---------------
1488  //-- NUCLEUS ---
1489  //---------------
1490  {int lz = 0 ;
1491  int nr = mgrid.get_nr(lz) ;
1492  double alpha = mp_aff->get_alpha()[lz] ;
1493  Matrice ope(2*nr, 2*nr) ;
1494  if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
1495 
1496  for (int k=0 ; k<np+1 ; k++) {
1497  for (int j=0 ; j<nt ; j++) {
1498  // quantic numbers and spectral bases
1499  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
1500  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q < 2)) {
1501  if (need_calculation) {
1502  ope.set_etat_qcq() ;
1503  Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
1504  Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
1505 
1506  for (int lin=0; lin<nr; lin++)
1507  for (int col=0; col<nr; col++)
1508  ope.set(lin,col) = md(lin,col) + 3*ms(lin,col) ;
1509  for (int lin=0; lin<nr; lin++)
1510  for (int col=0; col<nr; col++)
1511  ope.set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
1512  for (int lin=0; lin<nr; lin++)
1513  for (int col=0; col<nr; col++)
1514  ope.set(lin+nr,col) = -0.5*ms(lin,col) ;
1515  for (int lin=0; lin<nr; lin++)
1516  for (int col=0; col<nr; col++)
1517  ope.set(lin+nr,col+nr) = md(lin,col) + 3*ms(lin, col);
1518 
1519  ope *= 1./alpha ;
1520  for (int col=0; col<2*nr; col++) {
1521  ope.set(nr-1, col) = 0 ;
1522  ope.set(2*nr-1, col) = 0 ;
1523  }
1524  int pari = 1 ;
1525  if (base_r == R_CHEBP) {
1526  for (int col=0; col<nr; col++) {
1527  ope.set(nr-1, col) = pari ;
1528  ope.set(2*nr-1, col+nr) = pari ;
1529  pari = - pari ;
1530  }
1531  }
1532  else { //In the odd case, the last coefficient must be zero!
1533  ope.set(nr-1, nr-1) = 1 ;
1534  ope.set(2*nr-1, 2*nr-1) = 1 ;
1535  }
1536 
1537  ope.set_lu() ;
1538  if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
1539  Matrice* pope = new Matrice(ope) ;
1540  par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
1541  mat_done.set(l_q) = 1 ;
1542  }
1543  } //End of case when a calculation is needed
1544 
1545  const Matrice& oper = (par_mat == 0x0 ? ope :
1546  par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
1547  Tbl sec(2*nr) ;
1548  sec.set_etat_qcq() ;
1549  for (int lin=0; lin<nr; lin++)
1550  sec.set(lin) = (*source.get_spectral_va().c_cf)(lz, k, j, lin) ;
1551  for (int lin=0; lin<nr; lin++)
1552  sec.set(nr+lin) = -0.5*(*source.get_spectral_va().c_cf)
1553  (lz, k, j, lin) ;
1554  sec.set(nr-1) = 0 ;
1555  if (base_r == R_CHEBP) {
1556  double h0 = 0 ; //In the l=0 case: 3*hrr(r=0) = h(r=0)
1557  int pari = 1 ;
1558  for (int col=0; col<nr; col++) {
1559  h0 += pari*
1560  (*source_coq.get_spectral_va().c_cf)(lz, k, j, col) ;
1561  pari = - pari ;
1562  }
1563  sec.set(nr-1) = h0 / 3. ;
1564  }
1565  sec.set(2*nr-1) = 0 ;
1566  Tbl sol = oper.inverse(sec) ;
1567  for (int i=0; i<nr; i++) {
1568  sol_part_hrr.set(lz, k, j, i) = sol(i) ;
1569  sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
1570  }
1571  sec.annule_hard() ;
1572  }
1573  }
1574  }
1575  }
1576 
1577  //-------------
1578  // -- Shells --
1579  //-------------
1580 
1581  for (int lz=1; lz<nz-1; lz++) {
1582  if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
1583  int nr = mgrid.get_nr(lz) ;
1584  int ind0 = 0 ;
1585  int ind1 = nr ;
1586  assert(mgrid.get_nt(lz) == nt) ;
1587  assert(mgrid.get_np(lz) == np) ;
1588  double alpha = mp_aff->get_alpha()[lz] ;
1589  double ech = mp_aff->get_beta()[lz] / alpha ;
1590  Matrice ope(2*nr, 2*nr) ;
1591 
1592  for (int k=0 ; k<np+1 ; k++) {
1593  for (int j=0 ; j<nt ; j++) {
1594  // quantic numbers and spectral bases
1595  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
1596  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q < 2)) {
1597  if (need_calculation) {
1598  ope.set_etat_qcq() ;
1599  Diff_xdsdx oxd(base_r, nr) ; const Matrice& mxd = oxd.get_matrice() ;
1600  Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
1601  Diff_id oid(base_r, nr) ; const Matrice& mid = oid.get_matrice() ;
1602 
1603  for (int lin=0; lin<nr; lin++)
1604  for (int col=0; col<nr; col++)
1605  ope.set(lin,col) = mxd(lin,col) + ech*md(lin,col)
1606  + 3*mid(lin,col) ;
1607  for (int lin=0; lin<nr; lin++)
1608  for (int col=0; col<nr; col++)
1609  ope.set(lin,col+nr) = -l_q*(l_q+1)*mid(lin,col) ;
1610  for (int lin=0; lin<nr; lin++)
1611  for (int col=0; col<nr; col++)
1612  ope.set(lin+nr,col) = -0.5*mid(lin,col) ;
1613  for (int lin=0; lin<nr; lin++)
1614  for (int col=0; col<nr; col++)
1615  ope.set(lin+nr,col+nr) = mxd(lin,col) + ech*md(lin,col)
1616  + 3*mid(lin, col) ;
1617 
1618  for (int col=0; col<2*nr; col++) {
1619  ope.set(ind0+nr-1, col) = 0 ;
1620  ope.set(ind1+nr-1, col) = 0 ;
1621  }
1622  ope.set(ind0+nr-1, ind0) = 1 ;
1623  ope.set(ind1+nr-1, ind1) = 1 ;
1624 
1625  ope.set_lu() ;
1626  if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
1627  Matrice* pope = new Matrice(ope) ;
1628  par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
1629  mat_done.set(l_q) = 1 ;
1630  }
1631  } //End of case when a calculation is needed
1632  const Matrice& oper = (par_mat == 0x0 ? ope :
1633  par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
1634  Tbl sec(2*nr) ;
1635  sec.set_etat_qcq() ;
1636  for (int lin=0; lin<nr; lin++)
1637  sec.set(lin) = (*source_coq.get_spectral_va().c_cf)
1638  (lz, k, j, lin) ;
1639  for (int lin=0; lin<nr; lin++)
1640  sec.set(nr+lin) = -0.5*(*source_coq.get_spectral_va().c_cf)
1641  (lz, k, j, lin) ;
1642  sec.set(ind0+nr-1) = 0 ;
1643  sec.set(ind1+nr-1) = 0 ;
1644  Tbl sol = oper.inverse(sec) ;
1645 
1646  for (int i=0; i<nr; i++) {
1647  sol_part_hrr.set(lz, k, j, i) = sol(i) ;
1648  sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
1649  }
1650  sec.annule_hard() ;
1651  sec.set(ind0+nr-1) = 1 ;
1652  sol = oper.inverse(sec) ;
1653  for (int i=0; i<nr; i++) {
1654  sol_hom1_hrr.set(lz, k, j, i) = sol(i) ;
1655  sol_hom1_eta.set(lz, k, j, i) = sol(i+nr) ;
1656  }
1657  sec.set(ind0+nr-1) = 0 ;
1658  sec.set(ind1+nr-1) = 1 ;
1659  sol = oper.inverse(sec) ;
1660  for (int i=0; i<nr; i++) {
1661  sol_hom2_hrr.set(lz, k, j, i) = sol(i) ;
1662  sol_hom2_eta.set(lz, k, j, i) = sol(i+nr) ;
1663  }
1664  }
1665  }
1666  }
1667  }
1668 
1669  //------------------------------
1670  // Compactified external domain
1671  //------------------------------
1672  {int lz = nz-1 ;
1673  if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
1674  int nr = mgrid.get_nr(lz) ;
1675  int ind0 = 0 ;
1676  int ind1 = nr ;
1677  assert(mgrid.get_nt(lz) == nt) ;
1678  assert(mgrid.get_np(lz) == np) ;
1679  double alpha = mp_aff->get_alpha()[lz] ;
1680  Matrice ope(2*nr, 2*nr) ;
1681 
1682  for (int k=0 ; k<np+1 ; k++) {
1683  for (int j=0 ; j<nt ; j++) {
1684  // quantic numbers and spectral bases
1685  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
1686  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q < 2)) {
1687  if (need_calculation) {
1688  ope.set_etat_qcq() ;
1689  Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
1690  Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
1691 
1692  for (int lin=0; lin<nr; lin++)
1693  for (int col=0; col<nr; col++)
1694  ope.set(lin,col) = - md(lin,col) + 3*ms(lin,col) ;
1695  for (int lin=0; lin<nr; lin++)
1696  for (int col=0; col<nr; col++)
1697  ope.set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
1698  for (int lin=0; lin<nr; lin++)
1699  for (int col=0; col<nr; col++)
1700  ope.set(lin+nr,col) = -0.5*ms(lin,col) ;
1701  for (int lin=0; lin<nr; lin++)
1702  for (int col=0; col<nr; col++)
1703  ope.set(lin+nr,col+nr) = -md(lin,col) + 3*ms(lin, col) ;
1704 
1705  ope *= 1./alpha ;
1706  for (int col=0; col<2*nr; col++) {
1707  ope.set(ind0+nr-2, col) = 0 ;
1708  ope.set(ind0+nr-1, col) = 0 ;
1709  ope.set(ind1+nr-2, col) = 0 ;
1710  ope.set(ind1+nr-1, col) = 0 ;
1711  }
1712  for (int col=0; col<nr; col++) {
1713  ope.set(ind0+nr-1, col+ind0) = 1 ;
1714  ope.set(ind1+nr-1, col+ind1) = 1 ;
1715  }
1716  ope.set(ind0+nr-2, ind0+1) = 1 ;
1717  ope.set(ind1+nr-2, ind1+1) = 1 ;
1718 
1719  ope.set_lu() ;
1720  if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
1721  Matrice* pope = new Matrice(ope) ;
1722  par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
1723  mat_done.set(l_q) = 1 ;
1724  }
1725  } //End of case when a calculation is needed
1726  const Matrice& oper = (par_mat == 0x0 ? ope :
1727  par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
1728  Tbl sec(2*nr) ;
1729  sec.set_etat_qcq() ;
1730  for (int lin=0; lin<nr; lin++)
1731  sec.set(lin) = (*source.get_spectral_va().c_cf)
1732  (lz, k, j, lin) ;
1733  for (int lin=0; lin<nr; lin++)
1734  sec.set(nr+lin) = -0.5*(*source.get_spectral_va().c_cf)
1735  (lz, k, j, lin) ;
1736  sec.set(ind0+nr-2) = 0 ;
1737  sec.set(ind0+nr-1) = 0 ;
1738  sec.set(ind1+nr-2) = 0 ;
1739  sec.set(ind1+nr-1) = 0 ;
1740  Tbl sol = oper.inverse(sec) ;
1741  for (int i=0; i<nr; i++) {
1742  sol_part_hrr.set(lz, k, j, i) = sol(i) ;
1743  sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
1744  }
1745  sec.annule_hard() ;
1746  sec.set(ind0+nr-2) = 1 ;
1747  sol = oper.inverse(sec) ;
1748  for (int i=0; i<nr; i++) {
1749  sol_hom1_hrr.set(lz, k, j, i) = sol(i) ;
1750  sol_hom1_eta.set(lz, k, j, i) = sol(i+nr) ;
1751  }
1752  sec.set(ind0+nr-2) = 0 ;
1753  sec.set(ind1+nr-2) = 1 ;
1754  sol = oper.inverse(sec) ;
1755  for (int i=0; i<nr; i++) {
1756  sol_hom2_hrr.set(lz, k, j, i) = sol(i) ;
1757  sol_hom2_eta.set(lz, k, j, i) = sol(i+nr) ;
1758  }
1759  }
1760  }
1761  }
1762  }
1763 
1764  int taille = 2*(nz-1) ;
1765  Mtbl_cf& mhrr = *hrr.set_spectral_va().c_cf ;
1766  Mtbl_cf& meta = *tilde_eta.set_spectral_va().c_cf ;
1767 
1768  Tbl sec_membre(taille) ;
1769  Matrice systeme(taille, taille) ;
1770  int ligne ; int colonne ;
1771 
1772  // Loop on l and m
1773  //----------------
1774  for (int k=0 ; k<np+1 ; k++)
1775  for (int j=0 ; j<nt ; j++) {
1776  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
1777  if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q < 2)) {
1778  ligne = 0 ;
1779  colonne = 0 ;
1780  systeme.annule_hard() ;
1781  sec_membre.annule_hard() ;
1782 
1783  //Nucleus
1784  int nr = mgrid.get_nr(0) ;
1785 
1786  sec_membre.set(ligne) = -sol_part_hrr.val_out_bound_jk(0, j, k) ;
1787  ligne++ ;
1788 
1789  sec_membre.set(ligne) = -sol_part_eta.val_out_bound_jk(0, j, k) ;
1790 
1791  //shells
1792  for (int zone=1 ; zone<nz-1 ; zone++) {
1793  nr = mgrid.get_nr(zone) ;
1794  ligne-- ;
1795 
1796  //Condition at x = -1
1797  systeme.set(ligne, colonne) =
1798  - sol_hom1_hrr.val_in_bound_jk(zone, j, k) ;
1799  systeme.set(ligne, colonne+1) =
1800  - sol_hom2_hrr.val_in_bound_jk(zone, j, k) ;
1801 
1802  sec_membre.set(ligne) += sol_part_hrr.val_in_bound_jk(zone, j, k) ;
1803  ligne++ ;
1804 
1805  systeme.set(ligne, colonne) =
1806  - sol_hom1_eta.val_in_bound_jk(zone, j, k) ;
1807  systeme.set(ligne, colonne+1) =
1808  - sol_hom2_eta.val_in_bound_jk(zone, j, k) ;
1809 
1810  sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(zone, j, k) ;
1811  ligne++ ;
1812 
1813  // Condition at x=1
1814  systeme.set(ligne, colonne) =
1815  sol_hom1_hrr.val_out_bound_jk(zone, j, k) ;
1816  systeme.set(ligne, colonne+1) =
1817  sol_hom2_hrr.val_out_bound_jk(zone, j, k) ;
1818 
1819  sec_membre.set(ligne) -= sol_part_hrr.val_out_bound_jk(zone, j, k) ;
1820  ligne++ ;
1821 
1822  systeme.set(ligne, colonne) =
1823  sol_hom1_eta.val_out_bound_jk(zone, j, k) ;
1824  systeme.set(ligne, colonne+1) =
1825  sol_hom2_eta.val_out_bound_jk(zone, j, k) ;
1826 
1827  sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone, j, k) ;
1828 
1829  colonne += 2 ;
1830  }
1831 
1832  //Compactified external domain
1833  nr = mgrid.get_nr(nz-1) ;
1834 
1835  ligne-- ;
1836 
1837  systeme.set(ligne, colonne) =
1838  - sol_hom1_hrr.val_in_bound_jk(nz-1, j, k) ;
1839  systeme.set(ligne, colonne+1) =
1840  - sol_hom2_hrr.val_in_bound_jk(nz-1, j, k) ;
1841 
1842  sec_membre.set(ligne) += sol_part_hrr.val_in_bound_jk(nz-1, j, k) ;
1843  ligne++ ;
1844 
1845  systeme.set(ligne, colonne) =
1846  - sol_hom1_eta.val_in_bound_jk(nz-1, j, k) ;
1847  systeme.set(ligne, colonne+1) =
1848  - sol_hom2_eta.val_in_bound_jk(nz-1, j, k) ;
1849 
1850  sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(nz-1, j, k) ;
1851 
1852  // Solution of the system giving the coefficients for the homogeneous
1853  // solutions
1854  //-------------------------------------------------------------------
1855  systeme.set_lu() ;
1856  Tbl facteur = systeme.inverse(sec_membre) ;
1857  int conte = 0 ;
1858 
1859  // everything is put to the right place...
1860  //----------------------------------------
1861  nr = mgrid.get_nr(0) ; //nucleus
1862  for (int i=0 ; i<nr ; i++) {
1863  mhrr.set(0, k, j, i) = sol_part_hrr(0, k, j, i) ;
1864  meta.set(0, k, j, i) = sol_part_eta(0, k, j, i) ;
1865  }
1866  for (int zone=1 ; zone<nz-1 ; zone++) { //shells
1867  nr = mgrid.get_nr(zone) ;
1868  for (int i=0 ; i<nr ; i++) {
1869  mhrr.set(zone, k, j, i) = sol_part_hrr(zone, k, j, i)
1870  + facteur(conte)*sol_hom1_hrr(zone, k, j, i)
1871  + facteur(conte+1)*sol_hom2_hrr(zone, k, j, i) ;
1872 
1873  meta.set(zone, k, j, i) = sol_part_eta(zone, k, j, i)
1874  + facteur(conte)*sol_hom1_eta(zone, k, j, i)
1875  + facteur(conte+1)*sol_hom2_eta(zone, k, j, i) ;
1876  }
1877  conte+=2 ;
1878  }
1879  nr = mgrid.get_nr(nz-1) ; //compactified external domain
1880  for (int i=0 ; i<nr ; i++) {
1881  mhrr.set(nz-1, k, j, i) = sol_part_hrr(nz-1, k, j, i)
1882  + facteur(conte)*sol_hom1_hrr(nz-1, k, j, i)
1883  + facteur(conte+1)*sol_hom2_hrr(nz-1, k, j, i) ;
1884 
1885  meta.set(nz-1, k, j, i) = sol_part_eta(nz-1, k, j, i)
1886  + facteur(conte)*sol_hom1_eta(nz-1, k, j, i)
1887  + facteur(conte+1)*sol_hom2_eta(nz-1, k, j, i) ;
1888  }
1889 
1890  } // End of nullite_plm
1891  } //End of loop on theta
1892 }
1893 
1894 }
Bases of the spectral expansions.
Definition: base_val.h:322
void mult_x()
The basis is transformed as with a multiplication by .
int give_lmax(const Mg3d &mgrid, int lz) const
Returns the highest multipole for a given grid.
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:129
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_dsdx.C:94
Class for the elementary differential operator Identity (see the base class Diff ).
Definition: diff.h:210
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_id.C:91
Class for the elementary differential operator division by (see the base class Diff ).
Definition: diff.h:329
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_sx.C:100
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:409
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_xdsdx.C:98
Basic integer array class.
Definition: itbl.h:122
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: itbl.C:261
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
void annule_hard()
Sets the Itbl to zero in a hard way.
Definition: itbl.C:272
Affine radial mapping.
Definition: map.h:2027
const double * get_beta() const
Returns the pointer on the array beta.
Definition: map_af.C:481
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:477
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Matrix handling.
Definition: matrice.h:152
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:175
double & set(int j, int i)
Read/write of a particuliar element.
Definition: matrice.h:277
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition: matrice.C:424
void annule_hard()
Sets the logical state to ETATQCQ (undefined state).
Definition: matrice.C:193
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition: matrice.C:392
Multi-domain grid.
Definition: grilles.h:273
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Definition: grilles.h:485
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Definition: grilles.h:495
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:474
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:186
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition: mtbl_cf.C:312
double val_out_bound_jk(int l, int j, int k) const
Computes the angular coefficient of index j,k of the field represented by *this at by means of the s...
double val_in_bound_jk(int l, int j, int k) const
Computes the angular coefficient of index j,k of the field represented by *this at by means of the s...
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:294
Parameter storage.
Definition: param.h:125
int get_n_itbl_mod() const
Returns the number of modifiable Itbl 's addresses in the list.
Definition: param.C:722
int get_n_matrice_mod() const
Returns the number of modifiable Matrice 's addresses in the list.
Definition: param.C:859
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:292
int get_n_tbl_mod() const
Returns the number of modifiable Tbl 's addresses in the list.
Definition: param.C:584
Itbl & get_itbl_mod(int position=0) const
Returns the reference of a stored modifiable Itbl .
Definition: param.C:774
int get_n_int_mod() const
Returns the number of modifiable int 's addresses in the list.
Definition: param.C:378
Tbl & get_tbl_mod(int position=0) const
Returns the reference of a modifiable Tbl stored in the list.
Definition: param.C:636
void clean_all()
Deletes all the objects stored as modifiables, i.e.
Definition: param.C:174
void add_matrice_mod(Matrice &ti, int position=0)
Adds the address of a new modifiable Matrice to the list.
Definition: param.C:866
Matrice & get_matrice_mod(int position=0) const
Returns the reference of a modifiable Matrice stored in the list.
Definition: param.C:911
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 add_tbl_mod(Tbl &ti, int position=0)
Adds the address of a new modifiable Tbl to the list.
Definition: param.C:591
void add_itbl_mod(Itbl &ti, int position=0)
Adds the address of a new modifiable Itbl to the list.
Definition: param.C:729
int get_n_int() const
Returns the number of stored int 's addresses.
Definition: param.C:239
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Definition: param.C:430
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
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
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va
Definition: scalar.h:1294
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 set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:353
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition: scalar.C:873
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:380
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:601
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
void mult_r_dzpuis(int ced_mult_r)
Multiplication by r everywhere but with the output flag dzpuis set to ced_mult_r .
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
Definition: scalar.C:797
void sol_Dirac_tilde_B(const Scalar &tilde_b, const Scalar &hh, Scalar &hrr, Scalar &tilde_eta, Scalar &www, Param *par_bc=0x0, Param *par_mat=0x0) const
Solves a system of three coupled first-order PDEs obtained from divergence-free conditions (Dirac gau...
void sol_Dirac_A(const Scalar &aaa, Scalar &tilde_mu, Scalar &xxx, const Param *par_bc=0x0) const
Solves a system of two coupled first-order PDEs obtained from the divergence-free condition (Dirac ga...
void sol_Dirac_l01(const Scalar &hh, Scalar &hrr, Scalar &tilde_eta, Param *par_mat) const
Solves the same system as Sym_tensor_trans::sol_Dirac_tilde_B but only for .
Basic array class.
Definition: tbl.h:161
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:372
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
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
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
Mtbl * c
Values of the function at the points of the multi-grid
Definition: valeur.h:299
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:295
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:666
Lorene prototypes.
Definition: app_hor.h:64