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