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