28 char sol_Dirac_A_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Tensor/vector_divfree_A.C,v 1.6 2014/10/13 08:53:44 j_novak Exp $" ;
78 const Param* par_bc)
const {
81 assert(mp_aff != 0x0) ;
90 assert(aaa.
get_etat() != ETATNONDEF) ;
93 bool ced = (mgrid.
get_type_r(nzm1) == UNSURR) ;
94 int n_shell = ced ? nz-2 : nzm1 ;
98 n_shell = (nz_bc < n_shell ? nz_bc : n_shell) ;
99 bool cedbc = (ced && (nz_bc == nzm1)) ;
102 assert(par_bc != 0x0) ;
106 int nt = mgrid.
get_nt(0) ;
107 int np = mgrid.
get_np(0) ;
135 int l_q, m_q, base_r ;
141 int nr = mgrid.
get_nr(lz) ;
146 for (
int k=0 ; k<np+1 ; k++) {
147 for (
int j=0 ; j<nt ; j++) {
150 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 0)) {
154 for (
int lin=0; lin<nr; lin++)
155 for (
int col=0; col<nr; col++)
156 ope.
set(lin,col) = md(lin,col) + 2*ms(lin,col) ;
157 for (
int lin=0; lin<nr; lin++)
158 for (
int col=0; col<nr; col++)
159 ope.
set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
160 for (
int lin=0; lin<nr; lin++)
161 for (
int col=0; col<nr; col++)
162 ope.
set(lin+nr,col) = -ms(lin,col) ;
163 for (
int lin=0; lin<nr; lin++)
164 for (
int col=0; col<nr; col++)
165 ope.
set(lin+nr,col+nr) = md(lin,col) + ms(lin,col) ;
173 for (
int col=0 ; col<nr; col++) {
174 ope.
set(nr-1,col) = pari ;
175 ope.
set(nr-1,col+nr) = -pari ;
178 ope.
set(2*nr-1, nr) = 1;
181 for (
int col=0; col<2*nr; col++)
182 ope.
set(ind1+nr-2, col) = 0 ;
183 for (
int col=0; col<2*nr; col++) {
184 ope.
set(nr-1, col) = 0 ;
185 ope.
set(2*nr-1, col) = 0 ;
189 for (
int col=0; col<nr; col++) {
190 ope.
set(nr-1, col) = pari ;
191 ope.
set(2*nr-1, col+nr) = pari ;
196 ope.
set(nr-1, nr-1) = 1 ;
197 ope.
set(2*nr-1, 2*nr-1) = 1 ;
200 ope.
set(ind1+nr-2, ind1) = 1 ;
207 for (
int lin=0; lin<nr; lin++)
209 for (
int lin=0; lin<nr; lin++)
212 sec.
set(2*nr-1) = 0 ;
213 sec.
set(ind1+nr-2) = 0 ;
218 for (
int i=0; i<nr; i++) {
219 sol_part_vr.
set(lz, k, j, i) = sol(i) ;
220 sol_part_eta.
set(lz, k, j, i) = sol(i+nr) ;
223 sec.
set(ind1+nr-2) = 1 ;
227 for (
int i=0; i<nr; i++) {
228 sol_hom2_vr.
set(lz, k, j, i) = sol(i) ;
229 sol_hom2_eta.
set(lz, k, j, i) = sol(i+nr) ;
240 for (
int lz=1; lz <= n_shell; lz++) {
241 int nr = mgrid.
get_nr(lz) ;
242 assert(mgrid.
get_nt(lz) == nt) ;
243 assert(mgrid.
get_np(lz) == np) ;
245 double ech = mp_aff->
get_beta()[lz] / alpha ;
249 for (
int k=0 ; k<np+1 ; k++) {
250 for (
int j=0 ; j<nt ; j++) {
253 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 0)) {
258 for (
int lin=0; lin<nr; lin++)
259 for (
int col=0; col<nr; col++)
260 ope.
set(lin,col) = mxd(lin,col) + ech*md(lin,col)
262 for (
int lin=0; lin<nr; lin++)
263 for (
int col=0; col<nr; col++)
264 ope.
set(lin,col+nr) = -l_q*(l_q+1)*mid(lin,col) ;
265 for (
int lin=0; lin<nr; lin++)
266 for (
int col=0; col<nr; col++)
267 ope.
set(lin+nr,col) = -mid(lin,col) ;
268 for (
int lin=0; lin<nr; lin++)
269 for (
int col=0; col<nr; col++)
270 ope.
set(lin+nr,col+nr) = mxd(lin,col) + ech*md(lin,col) + mid(lin,col) ;
274 for (
int col=0; col<2*nr; col++) {
275 ope.
set(ind0+nr-1, col) = 0 ;
276 ope.
set(ind1+nr-1, col) = 0 ;
278 ope.
set(ind0+nr-1, ind0) = 1 ;
279 ope.
set(ind1+nr-1, ind1) = 1 ;
285 for (
int lin=0; lin<nr; lin++)
287 for (
int lin=0; lin<nr; lin++)
290 sec.
set(ind0+nr-1) = 0 ;
291 sec.
set(ind1+nr-1) = 0 ;
293 for (
int i=0; i<nr; i++) {
294 sol_part_vr.
set(lz, k, j, i) = sol(i) ;
295 sol_part_eta.
set(lz, k, j, i) = sol(i+nr) ;
300 sec.
set(ind0+nr-1) = 1 ;
304 for (
int i=0; i<nr; i++) {
305 sol_hom1_vr.
set(lz, k, j, i) = sol(i) ;
306 sol_hom1_eta.
set(lz, k, j, i) = sol(i+nr) ;
308 sec.
set(ind0+nr-1) = 0 ;
309 sec.
set(ind1+nr-1) = 1 ;
314 for (
int i=0; i<nr; i++) {
315 sol_hom2_vr.
set(lz, k, j, i) = sol(i) ;
316 sol_hom2_eta.
set(lz, k, j, i) = sol(i+nr) ;
326 if (cedbc) {
int lz = nzm1 ;
327 int nr = mgrid.
get_nr(lz) ;
328 assert(mgrid.
get_nt(lz) == nt) ;
329 assert(mgrid.
get_np(lz) == np) ;
334 for (
int k=0 ; k<np+1 ; k++) {
335 for (
int j=0 ; j<nt ; j++) {
338 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 0)) {
343 for (
int lin=0; lin<nr; lin++)
344 for (
int col=0; col<nr; col++)
345 ope.
set(lin,col) = - md(lin,col) + 2*ms(lin,col) ;
346 for (
int lin=0; lin<nr; lin++)
347 for (
int col=0; col<nr; col++)
348 ope.
set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
349 for (
int lin=0; lin<nr; lin++)
350 for (
int col=0; col<nr; col++)
351 ope.
set(lin+nr,col) = -ms(lin,col) ;
352 for (
int lin=0; lin<nr; lin++)
353 for (
int col=0; col<nr; col++)
354 ope.
set(lin+nr,col+nr) = -md(lin,col) + ms(lin,col) ;
359 for (
int col=0; col<2*nr; col++) {
360 ope.
set(ind0+nr-1, col) = 0 ;
361 ope.
set(ind1+nr-2, col) = 0 ;
362 ope.
set(ind1+nr-1, col) = 0 ;
364 for (
int col=0; col<nr; col++) {
365 ope.
set(ind0+nr-1, col+ind0) = 1 ;
366 ope.
set(ind1+nr-1, col+ind1) = 1 ;
368 ope.
set(ind1+nr-2, ind1+1) = 1 ;
374 for (
int lin=0; lin<nr; lin++)
376 for (
int lin=0; lin<nr; lin++)
379 sec.
set(ind0+nr-1) = 0 ;
380 sec.
set(ind1+nr-2) = 0 ;
381 sec.
set(ind1+nr-1) = 0 ;
384 for (
int i=0; i<nr; i++) {
385 sol_part_vr.
set(lz, k, j, i) = sol(i) ;
386 sol_part_eta.
set(lz, k, j, i) = sol(i+nr) ;
389 sec.
set(ind1+nr-2) = 1 ;
391 for (
int i=0; i<nr; i++) {
392 sol_hom1_vr.
set(lz, k, j, i) = sol(i) ;
393 sol_hom1_eta.
set(lz, k, j, i) = sol(i+nr) ;
400 int taille = 2*nz_bc + 1;
401 if (cedbc) taille-- ;
405 Tbl sec_membre(taille) ;
406 Matrice systeme(taille, taille) ;
407 int ligne ;
int colonne ;
410 double c_vr = (cedbc ? 0 : par_bc->
get_tbl_mod(0)(0) ) ;
411 double d_vr = (cedbc ? 0 : par_bc->
get_tbl_mod(0)(1) ) ;
412 double c_eta = (cedbc ? 0 : par_bc->
get_tbl_mod(0)(2) ) ;
413 double d_eta = (cedbc ? 0 : par_bc->
get_tbl_mod(0)(3) ) ;
417 Mtbl_cf dhom1_vr = sol_hom1_vr ;
418 Mtbl_cf dhom2_vr = sol_hom2_vr ;
419 Mtbl_cf dpart_vr = sol_part_vr ;
420 Mtbl_cf dhom1_eta = sol_hom1_eta ;
421 Mtbl_cf dhom2_eta = sol_hom2_eta ;
422 Mtbl_cf dpart_eta = sol_part_eta ;
435 for (
int k=0 ; k<np+1 ; k++)
436 for (
int j=0 ; j<nt ; j++) {
438 if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q > 0)) {
447 double alpha = mp_aff->
get_alpha()[nz_bc] ;
448 systeme.
set(ligne, colonne) =
472 for (
int zone=1 ; zone<nz_bc ; zone++) {
477 systeme.
set(ligne, colonne) =
479 systeme.
set(ligne, colonne+1) =
485 systeme.
set(ligne, colonne) =
487 systeme.
set(ligne, colonne+1) =
494 systeme.
set(ligne, colonne) =
496 systeme.
set(ligne, colonne+1) =
502 systeme.
set(ligne, colonne) =
504 systeme.
set(ligne, colonne+1) =
513 nr = mgrid.
get_nr(nz_bc) ;
514 double alpha = mp_aff->
get_alpha()[nz_bc] ;
518 systeme.
set(ligne, colonne) =
520 if (!cedbc) systeme.
set(ligne, colonne+1) =
526 systeme.
set(ligne, colonne) =
528 if (!cedbc) systeme.
set(ligne, colonne+1) =
535 systeme.
set(ligne, colonne) =
540 systeme.
set(ligne, colonne+1) =
565 for (
int i=0 ; i<nr ; i++) {
566 mvr.
set(0, k, j, i) = sol_part_vr(0, k, j, i)
567 + facteur(conte)*sol_hom2_vr(0, k, j, i) ;
568 meta.
set(0, k, j, i) = sol_part_eta(0, k, j, i)
569 + facteur(conte)*sol_hom2_eta(0, k, j, i) ;
572 for (
int zone=1 ; zone<=n_shell ; zone++) {
574 for (
int i=0 ; i<nr ; i++) {
575 mvr.
set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
576 + facteur(conte)*sol_hom1_vr(zone, k, j, i)
577 + facteur(conte+1)*sol_hom2_vr(zone, k, j, i) ;
579 meta.
set(zone, k, j, i) = sol_part_eta(zone, k, j, i)
580 + facteur(conte)*sol_hom1_eta(zone, k, j, i)
581 + facteur(conte+1)*sol_hom2_eta(zone, k, j, i) ;
587 for (
int i=0 ; i<nr ; i++) {
588 mvr.
set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
589 + facteur(conte)*sol_hom1_vr(nzm1, k, j, i) ;
591 meta.
set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
592 + facteur(conte)*sol_hom1_eta(nzm1, k, j, i) ;