67#include "param_elliptic.h"
84 assert(mp_aff != 0x0) ;
109 assert(aaa.
get_etat() != ETATNONDEF) ;
112 bool ced = (mgrid.
get_type_r(nzm1) == UNSURR) ;
113 int n_shell = ced ? nz-2 : nzm1 ;
117 n_shell = (nz_bc < n_shell ? nz_bc : n_shell) ;
118 bool cedbc = (ced && (nz_bc == nzm1)) ;
121 assert(par_bc != 0x0) ;
125 int nt = mgrid.
get_nt(0) ;
126 int np = mgrid.
get_np(0) ;
154 int l_q, m_q, base_r ;
163 int nr = mgrid.
get_nr(lz) ;
168 for (
int k=0 ; k<np+1 ; k++) {
169 for (
int j=0 ; j<nt ; j++) {
172 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
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 ;
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 ;
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) ;
199 double ech = mp_aff->
get_beta()[lz] / alpha ;
203 for (
int k=0 ; k<np+1 ; k++) {
204 for (
int j=0 ; j<nt ; j++) {
207 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
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)
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) ;
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 ;
232 ope.
set(ind0+nr-1, ind0) = 1 ;
233 ope.
set(ind1+nr-1, ind1) = 1 ;
239 for (
int lin=0; lin<nr; lin++)
241 for (
int lin=0; lin<nr; lin++)
244 sec.
set(ind0+nr-1) = 0 ;
245 sec.
set(ind1+nr-1) = 0 ;
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) ;
252 sec.
set(ind0+nr-1) = 1 ;
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) ;
258 sec.
set(ind0+nr-1) = 0 ;
259 sec.
set(ind1+nr-1) = 1 ;
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) ;
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) ;
281 for (
int k=0 ; k<np+1 ; k++) {
282 for (
int j=0 ; j<nt ; j++) {
285 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
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) ;
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 ;
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 ;
314 ope.
set(ind1+nr-2, ind1+1) = 1 ;
320 for (
int lin=0; lin<nr; lin++)
322 for (
int lin=0; lin<nr; lin++)
325 sec.
set(ind0+nr-1) = 0 ;
326 sec.
set(ind1+nr-2) = 0 ;
327 sec.
set(ind1+nr-1) = 0 ;
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) ;
334 sec.
set(ind1+nr-2) = 1 ;
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) ;
348 int taille = 2*nz_bc ;
349 if (cedbc) taille-- ;
353 Tbl sec_membre(taille) ;
354 Matrice systeme(taille, taille) ;
355 int ligne ;
int colonne ;
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 ;
380 for (
int k=0 ; k<np+1 ; k++)
381 for (
int j=0 ; j<nt ; j++) {
383 if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
391 int nr = mgrid.
get_nr(1) ;
404 systeme.
set(ligne, colonne) =
406 systeme.
set(ligne, colonne+1) =
412 systeme.
set(ligne, colonne) =
414 systeme.
set(ligne, colonne+1) =
422 for (
int zone=2 ; zone<nz_bc ; zone++) {
427 systeme.
set(ligne, colonne) =
429 systeme.
set(ligne, colonne+1) =
435 systeme.
set(ligne, colonne) =
437 systeme.
set(ligne, colonne+1) =
444 systeme.
set(ligne, colonne) =
446 systeme.
set(ligne, colonne+1) =
452 systeme.
set(ligne, colonne) =
454 systeme.
set(ligne, colonne+1) =
463 nr = mgrid.
get_nr(nz_bc) ;
464 double alpha = mp_aff->
get_alpha()[nz_bc] ;
468 systeme.
set(ligne, colonne) =
470 if (!cedbc) systeme.
set(ligne, colonne+1) =
476 systeme.
set(ligne, colonne) =
478 if (!cedbc) systeme.
set(ligne, colonne+1) =
486 systeme.
set(ligne, colonne) =
491 systeme.
set(ligne, colonne+1) =
520 for (
int i=0 ; i<nr ; i++) {
521 mmu.
set(0, k, j, i) = 0 ;
522 mw.
set(0, k, j, i) = 0 ;
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);
534 for (
int zone=2 ; zone<=n_shell ; 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) ;
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) ;
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) ;
553 mw.
set(nzm1, k, j, i) = sol_part_x(nzm1, k, j, i)
554 + facteur(conte)*sol_hom1_x(nzm1, k, j, i) ;
584 assert(mp_aff != 0x0) ;
591 bool ced = (mgrid.
get_type_r(nzm1) == UNSURR) ;
592 int n_shell = ced ? nz-2 : nzm1 ;
597 n_shell = (nz_bc < n_shell ? nz_bc : n_shell) ;
598 bool cedbc = (ced && (nz_bc == nzm1)) ;
601 assert(par_bc != 0x0) ;
605 int nt = mgrid.
get_nt(0) ;
606 int np = mgrid.
get_np(0) ;
608 int l_q, m_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++) {
632 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1))
634 for (
int i=0; i<nr; i++)
657 Scalar tilde_b2 = tilde_b;
660 assert (tilde_b.
get_etat() != ETATNONDEF) ;
661 assert (hh.
get_etat() != ETATNONDEF) ;
664 Scalar source_coq = tilde_b ;
670 bool bnull = (tilde_b.
get_etat() == ETATZERO) ;
683 bool hnull = (hh.
get_etat() == ETATZERO) ;
690 bool need_calculation = true ;
691 if (par_mat != 0x0) {
692 bool param_new = false ;
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 ;
704 for (
int l=1; l<= n_shell; l++) {
707 mp_aff->
get_alpha()[l]) > 2.e-15) param_new = true ;
720 int* nz_bc_new =
new int(nz_bc) ;
722 int* lmax_new =
new int(lmax) ;
724 int* type_t_new =
new int(mgrid.
get_type_t()) ;
726 int* type_p_new =
new int(mgrid.
get_type_p()) ;
731 for (
int l=0; l<nz; l++)
733 Tbl* palpha =
new Tbl(nz) ;
737 for (
int l=1; l<nzm1; l++)
741 else need_calculation = false ;
773 Itbl mat_done(lmax) ;
783 int nr = mgrid.
get_nr(lz) ;
787 if (need_calculation && (par_mat != 0x0)) mat_done.
annule_hard() ;
789 for (
int k=0 ; k<np+1 ; k++) {
790 for (
int j=0 ; j<nt ; j++) {
793 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
794 if (need_calculation) {
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) ;
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 ;
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 ;
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 ;
852 ope.
set(ind2+nr-2, ind2) = 1 ;
855 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
858 mat_done.
set(l_q) = 1 ;
862 const Matrice& oper = (par_mat == 0x0 ? ope :
867 for (
int lin=0; lin<2*nr; lin++)
869 for (
int lin=0; lin<nr; lin++)
874 for (
int lin=0; lin<nr; lin++)
876 for (
int lin=0; lin<nr; lin++)
880 for (
int lin=0; lin<nr; lin++)
881 sec.
set(2*nr+lin) = -0.5/double(l_q+1)*(
886 for (
int lin=0; lin<nr; lin++)
887 sec.
set(2*nr+lin) = -0.5/double(l_q+1)*(
893 if (l_q>2) sec.
set(ind2+nr-2) = 0 ;
894 sec.
set(3*nr-1) = 0 ;
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) ;
903 sec.
set(ind2+nr-2) = 1 ;
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) ;
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) ;
934 double ech = mp_aff->
get_beta()[lz] / alpha ;
937 for (
int k=0 ; k<np+1 ; k++) {
938 for (
int j=0 ; j<nt ; j++) {
941 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
942 if (need_calculation) {
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)
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)
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 ;
986 ope.
set(ind0+nr-1, ind0) = 1 ;
987 ope.
set(ind1+nr-1, ind1) = 1 ;
988 ope.
set(ind2+nr-1, ind2) = 1 ;
991 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
994 mat_done.
set(l_q) = 1 ;
997 const Matrice& oper = (par_mat == 0x0 ? ope :
1002 for (
int lin=0; lin<2*nr; lin++)
1004 for (
int lin=0; lin<nr; lin++)
1009 for (
int lin=0; lin<nr; lin++)
1011 for (
int lin=0; lin<nr; lin++)
1015 for (
int lin=0; lin<nr; lin++)
1016 sec.
set(2*nr+lin) = -0.5/double(l_q+1)*(
1021 for (
int lin=0; lin<nr; lin++)
1022 sec.
set(2*nr+lin) = -0.5/double(l_q+1)*(
1028 sec.
set(ind0+nr-1) = 0 ;
1029 sec.
set(ind1+nr-1) = 0 ;
1030 sec.
set(ind2+nr-1) = 0 ;
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) ;
1038 sec.
set(ind0+nr-1) = 1 ;
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) ;
1045 sec.
set(ind0+nr-1) = 0 ;
1046 sec.
set(ind1+nr-1) = 1 ;
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) ;
1053 sec.
set(ind1+nr-1) = 0 ;
1054 sec.
set(ind2+nr-1) = 1 ;
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) ;
1069 if (cedbc) {
int lz = nzm1 ;
1070 if (need_calculation && (par_mat != 0x0)) mat_done.
annule_hard() ;
1071 int nr = mgrid.
get_nr(lz) ;
1075 double alpha = mp_aff->
get_alpha()[lz] ;
1078 for (
int k=0 ; k<np+1 ; k++) {
1079 for (
int j=0 ; j<nt ; j++) {
1082 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
1083 if (need_calculation) {
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) ;
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 ;
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 ;
1130 ope.
set(ind0+nr-2, ind0+1) = 1 ;
1131 ope.
set(ind1+nr-2, ind1+2) = 1 ;
1134 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
1137 mat_done.
set(l_q) = 1 ;
1140 const Matrice& oper = (par_mat == 0x0 ? ope :
1146 for (
int lin=0; lin<2*nr; lin++)
1148 for (
int lin=0; lin<nr; lin++)
1153 for (
int lin=0; lin<nr; lin++)
1155 for (
int lin=0; lin<nr; lin++)
1159 for (
int lin=0; lin<nr; lin++)
1160 sec.
set(2*nr+lin) = -0.5/double(l_q+1)*(
1165 for (
int lin=0; lin<nr; lin++)
1166 sec.
set(2*nr+lin) = -0.5/double(l_q+1)*(
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 ;
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) ;
1184 sec.
set(ind0+nr-2) = 1 ;
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) ;
1191 sec.
set(ind0+nr-2) = 0 ;
1192 sec.
set(ind1+nr-2) = 1 ;
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) ;
1208 int taille = 3*nz_bc ;
1209 if (cedbc) taille-- ;
1213 Tbl sec_membre(taille) ;
1214 Matrice systeme(taille, taille) ;
1215 int ligne ;
int colonne ;
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 ;
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();
1255 for (
int k=0 ; k<np+1 ; k++)
1256 for (
int j=0 ; j<nt ; j++) {
1258 if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
1266 int nr = mgrid.
get_nr(1);
1267 double alpha2 = mp_aff->
get_alpha()[1] ;
1270 const Coord& rr = (*mp_aff).r ;
1271 Scalar rrr(*mp_aff); rrr = rr;
1283 sec_membre.
set(ligne) += (*Bbb).val_in_bound_jk(1,j,k);
1290 systeme.
set(ligne, colonne) =
1292 systeme.
set(ligne, colonne+1) =
1294 systeme.
set(ligne, colonne+2) =
1298 double blob = (*(bound_eta.
get_spectral_va().c_cf)).val_in_bound_jk(1, j, k);
1306 systeme.
set(ligne, colonne) =
1308 systeme.
set(ligne, colonne+1) =
1310 systeme.
set(ligne, colonne+2) =
1316 systeme.
set(ligne, colonne) =
1318 systeme.
set(ligne, colonne+1) =
1320 systeme.
set(ligne, colonne+2) =
1326 systeme.
set(ligne, colonne) =
1328 systeme.
set(ligne, colonne+1) =
1330 systeme.
set(ligne, colonne+2) =
1338 for (
int zone=2 ; zone<nz_bc ; zone++) {
1339 nr = mgrid.
get_nr(zone) ;
1343 systeme.
set(ligne, colonne) =
1345 systeme.
set(ligne, colonne+1) =
1347 systeme.
set(ligne, colonne+2) =
1353 systeme.
set(ligne, colonne) =
1355 systeme.
set(ligne, colonne+1) =
1357 systeme.
set(ligne, colonne+2) =
1363 systeme.
set(ligne, colonne) =
1365 systeme.
set(ligne, colonne+1) =
1367 systeme.
set(ligne, colonne+2) =
1374 systeme.
set(ligne, colonne) =
1376 systeme.
set(ligne, colonne+1) =
1378 systeme.
set(ligne, colonne+2) =
1384 systeme.
set(ligne, colonne) =
1386 systeme.
set(ligne, colonne+1) =
1388 systeme.
set(ligne, colonne+2) =
1394 systeme.
set(ligne, colonne) =
1396 systeme.
set(ligne, colonne+1) =
1398 systeme.
set(ligne, colonne+2) =
1407 nr = mgrid.
get_nr(nz_bc) ;
1408 double alpha = mp_aff->
get_alpha()[nz_bc] ;
1412 systeme.
set(ligne, colonne) =
1414 systeme.
set(ligne, colonne+1) =
1416 if (!cedbc) systeme.
set(ligne, colonne+2) =
1422 systeme.
set(ligne, colonne) =
1424 systeme.
set(ligne, colonne+1) =
1426 if (!cedbc) systeme.
set(ligne, colonne+2) =
1432 systeme.
set(ligne, colonne) =
1434 systeme.
set(ligne, colonne+1) =
1436 if (!cedbc) systeme.
set(ligne, colonne+2) =
1443 systeme.
set(ligne, colonne) =
1450 systeme.
set(ligne, colonne+1) =
1457 systeme.
set(ligne, colonne+2) =
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 ;
1496 for (
int zone=1 ; zone<=n_shell ; zone++) {
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) ;
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) ;
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) ;
1517 nr = mgrid.
get_nr(nzm1) ;
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) ;
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) ;
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) ;
1559void Sym_tensor_trans::sol_Dirac_l01_2(
const Scalar& hh,
Scalar& hrr,
Scalar& tilde_eta,
1569 assert(mp_aff != 0x0) ;
1582 int nt = mgrid.
get_nt(0) ;
1583 int np = mgrid.
get_np(0) ;
1585 Scalar source = hh ;
1587 Scalar source_coq = hh ;
1589 source_coq.set_spectral_va().ylm() ;
1590 Base_val base = source.get_spectral_base() ;
1592 int lmax = base.give_lmax(mgrid, 0) + 1;
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() ;
1606 bool need_calculation = true ;
1611 int l_q, m_q, base_r ;
1612 Itbl mat_done(lmax) ;
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() ;
1623 for (
int k=0 ; k<np+1 ; k++) {
1624 for (
int j=0 ; j<nt ; j++) {
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() ;
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);
1647 for (
int col=0; col<2*nr; col++) {
1648 ope.set(nr-1, col) = 0 ;
1649 ope.set(2*nr-1, col) = 0 ;
1653 for (
int col=0; col<nr; col++) {
1654 ope.set(nr-1, col) = pari ;
1655 ope.set(2*nr-1, col+nr) = pari ;
1660 ope.set(nr-1, nr-1) = 1 ;
1661 ope.set(2*nr-1, 2*nr-1) = 1 ;
1665 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
1666 Matrice* pope =
new Matrice(ope) ;
1668 mat_done.set(l_q) = 1 ;
1672 const Matrice& oper = (par_mat == 0x0 ? ope :
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)
1685 for (
int col=0; col<nr; col++) {
1687 (*source_coq.get_spectral_va().c_cf)(lz, k, j, col) ;
1690 sec.set(nr-1) = h0 / 3. ;
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) ;
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) ;
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) ;
1719 for (
int k=0 ; k<np+1 ; k++) {
1720 for (
int j=0 ; j<nt ; j++) {
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() ;
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)
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)
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 ;
1749 ope.set(ind0+nr-1, ind0) = 1 ;
1750 ope.set(ind1+nr-1, ind1) = 1 ;
1753 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
1754 Matrice* pope =
new Matrice(ope) ;
1756 mat_done.set(l_q) = 1 ;
1759 const Matrice& oper = (par_mat == 0x0 ? ope :
1762 sec.set_etat_qcq() ;
1763 for (
int lin=0; lin<nr; lin++)
1764 sec.set(lin) = (*source_coq.get_spectral_va().c_cf)
1766 for (
int lin=0; lin<nr; lin++)
1767 sec.set(nr+lin) = -0.5*(*source_coq.get_spectral_va().c_cf)
1769 sec.set(ind0+nr-1) = 0 ;
1770 sec.set(ind1+nr-1) = 0 ;
1771 Tbl sol = oper.inverse(sec) ;
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) ;
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) ;
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) ;
1800 if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
1801 int nr = mgrid.
get_nr(lz) ;
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) ;
1809 for (
int k=0 ; k<np+1 ; k++) {
1810 for (
int j=0 ; j<nt ; j++) {
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() ;
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) ;
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 ;
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 ;
1843 ope.set(ind0+nr-2, ind0+1) = 1 ;
1844 ope.set(ind1+nr-2, ind1+1) = 1 ;
1847 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
1848 Matrice* pope =
new Matrice(ope) ;
1850 mat_done.set(l_q) = 1 ;
1853 const Matrice& oper = (par_mat == 0x0 ? ope :
1856 sec.set_etat_qcq() ;
1857 for (
int lin=0; lin<nr; lin++)
1858 sec.set(lin) = (*source.get_spectral_va().c_cf)
1860 for (
int lin=0; lin<nr; lin++)
1861 sec.set(nr+lin) = -0.5*(*source.get_spectral_va().c_cf)
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) ;
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) ;
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) ;
1891 int taille = 2*(nz-1) ;
1895 Tbl sec_membre(taille) ;
1896 Matrice systeme(taille, taille) ;
1897 int ligne ;
int colonne ;
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)) {
1907 systeme.annule_hard() ;
1908 sec_membre.annule_hard() ;
1911 int nr = mgrid.
get_nr(0) ;
1913 sec_membre.set(ligne) = -sol_part_hrr.val_out_bound_jk(0, j, k) ;
1916 sec_membre.set(ligne) = -sol_part_eta.val_out_bound_jk(0, j, k) ;
1919 for (
int zone=1 ; zone<nz-1 ; zone++) {
1920 nr = mgrid.
get_nr(zone) ;
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) ;
1929 sec_membre.set(ligne) += sol_part_hrr.val_in_bound_jk(zone, j, k) ;
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) ;
1937 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(zone, j, k) ;
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) ;
1946 sec_membre.set(ligne) -= sol_part_hrr.val_out_bound_jk(zone, j, k) ;
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) ;
1954 sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone, j, k) ;
1960 nr = mgrid.
get_nr(nz-1) ;
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) ;
1969 sec_membre.set(ligne) += sol_part_hrr.val_in_bound_jk(nz-1, j, k) ;
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) ;
1977 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(nz-1, j, k) ;
1983 Tbl facteur = systeme.inverse(sec_membre) ;
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) ;
1993 for (
int zone=1 ; zone<nz-1 ; zone++) {
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) ;
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) ;
2006 nr = mgrid.
get_nr(nz-1) ;
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) ;
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) ;
Bases of the spectral expansions.
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.
Class for the elementary differential operator (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Class for the elementary differential operator Identity (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Class for the elementary differential operator division by (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Class for the elementary differential operator (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Basic integer array class.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
int & set(int i)
Read/write of a particular element (index i ) (1D case)
void annule_hard()
Sets the Itbl to zero in a hard way.
const double * get_beta() const
Returns the pointer on the array beta.
const double * get_alpha() const
Returns the pointer on the array alpha.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double & set(int j, int i)
Read/write of a particuliar element.
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
void annule_hard()
Sets the logical state to ETATQCQ (undefined state).
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
int get_nzone() const
Returns the number of domains.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Coefficients storage for the multi-domain spectral method.
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
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...
int get_n_itbl_mod() const
Returns the number of modifiable Itbl 's addresses in the list.
int get_n_matrice_mod() const
Returns the number of modifiable Matrice 's addresses in the list.
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
int get_n_tbl_mod() const
Returns the number of modifiable Tbl 's addresses in the list.
Itbl & get_itbl_mod(int position=0) const
Returns the reference of a stored modifiable Itbl .
int get_n_int_mod() const
Returns the number of modifiable int 's addresses in the list.
Tbl & get_tbl_mod(int position=0) const
Returns the reference of a modifiable Tbl stored in the list.
void clean_all()
Deletes all the objects stored as modifiables, i.e.
void add_matrice_mod(Matrice &ti, int position=0)
Adds the address of a new modifiable Matrice to the list.
Matrice & get_matrice_mod(int position=0) const
Returns the reference of a modifiable Matrice stored in the list.
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
void add_tbl_mod(Tbl &ti, int position=0)
Adds the address of a new modifiable Tbl to the list.
void add_itbl_mod(Itbl &ti, int position=0)
Adds the address of a new modifiable Itbl to the list.
int get_n_int() const
Returns the number of stored int 's addresses.
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Tensor field of valence 0 (or component of a tensorial field).
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.
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).
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
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...
const Scalar & dsdr() const
Returns of *this .
Valeur & set_spectral_va()
Returns va (read/write version)
const Valeur & get_spectral_va() const
Returns va (read only version)
void annule_hard()
Sets the Scalar to zero in a hard way.
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va
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
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 .
void annule_hard()
Sets the Tbl to zero in a hard way.
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double & set(int i)
Read/write of a particular element (index i) (1D case)
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
void ylm()
Computes the coefficients of *this.
Mtbl * c
Values of the function at the points of the multi-grid
void coef_i() const
Computes the physical value of *this.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
void ylm_i()
Inverse of ylm()
#define R_CHEBP
base de Cheb. paire (rare) seulement
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
void annule_domain(int l)
Sets the Tensor to zero in a given domain.