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