23 char vector_divfree_A_tau[] =
"$Header: /cvsroot/Lorene/C++/Source/Tensor/vector_divfree_A_tau.C,v 1.4 2014/10/13 08:53:45 j_novak Exp $" ;
69 #include "type_parite.h"
76 assert(mp_aff != 0x0) ;
85 assert(aaa.
get_etat() != ETATNONDEF) ;
88 bool ced = (mgrid.
get_type_r(nzm1) == UNSURR) ;
89 int n_shell = ced ? nz-2 : nzm1 ;
93 n_shell = (nz_bc < n_shell ? nz_bc : n_shell) ;
102 int nt = mgrid.
get_nt(0) ;
103 int np = mgrid.
get_np(0) ;
137 for (
int l=0 ; l<nz ; l++) {
146 Tbl sec_membre (size) ;
153 int l_quant, m_quant ;
156 for (
int k=0 ; k<np+1 ; k++)
157 for (
int j=0 ; j<nt ; j++)
158 if (nullite_plm(j, nt, k, np, base) == 1) {
163 int column_courant = 0 ;
164 int ligne_courant = 0 ;
171 alpha = (*mp_aff).get_alpha()[0] ;
186 for (
int col=0 ; col<nr ; col++) {
188 systeme.
set(ligne_courant, col+column_courant) = 1 ;
189 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ; }
191 systeme.
set(ligne_courant, col+column_courant) = -1 ;
192 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ; }
197 for (
int col=0 ; col<nr ; col++) {
199 systeme.
set(ligne_courant, col+column_courant) = 2*col+1 ;
200 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ; }
202 systeme.
set(ligne_courant, col+column_courant) = -(2*col+1) ;
203 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ; }
215 for (
int col=0 ; col<nr ; col++) {
216 systeme.
set(ligne_courant, col+column_courant) = col*(col+1)*(col+2)*(col+3)/
double(12)*(2*(col%2)-1);
217 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ;
220 else if (l_quant == 1) {
222 for (
int col=0 ; col<nr ; col++) {
223 systeme.
set(ligne_courant, col+column_courant) = (col+1)*(col+2)/double(2)*(1-2*(col%2)) ;
224 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ;
229 for (
int col=0 ; col<nr ; col++) {
230 systeme.
set(ligne_courant, col+column_courant) = (col+1)*(col+2)/double(2)*(1-2*(col%2)) ;
231 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ;
232 systeme.
set(ligne_courant+1, col+column_courant) = col*(col+1)*(col+2)*(col+3)/
double(12)*(2*(col%2)-1) ;
233 systeme.
set(ligne_courant+1, col+column_courant+nr) = 0 ;
237 ligne_courant += nbr_cl ;
240 for (
int lig=0 ; lig<nr-1-nbr_cl ; lig++) {
241 for (
int col=0 ; col<nr ; col++) {
242 systeme.
set(lig+ligne_courant,col+column_courant) = mdn(lig,col) + msn(lig,col) ;
244 systeme.
set(lig+ligne_courant,col+column_courant+nr)= - msn(lig,col) ;
250 ligne_courant += nr-1-nbr_cl ;
260 for (
int col=0 ; col<nr ; col++) {
262 systeme.
set(ligne_courant, col+column_courant) = 0 ;
263 systeme.
set(ligne_courant, col+column_courant+nr) = 1 ; }
265 systeme.
set(ligne_courant, col+column_courant) = 0 ;
266 systeme.
set(ligne_courant, col+column_courant+nr) = -1 ; }
272 for (
int col=0 ; col<nr ; col++) {
274 systeme.
set(ligne_courant, col+column_courant) = 0 ;
275 systeme.
set(ligne_courant, col+column_courant+nr) = 2*col+1 ; }
277 systeme.
set(ligne_courant, col+column_courant) = 0 ;
278 systeme.
set(ligne_courant, col+column_courant+nr) = -(2*col+1) ; }
290 for (
int col=0 ; col<nr ; col++) {
291 systeme.
set(ligne_courant, col+column_courant) = 0 ;
292 systeme.
set(ligne_courant, col+column_courant+nr) = col*(col+1)*(col+2)*(col+3)/
double(12)*(2*(col%2)-1) ;
295 else if (l_quant == 1) {
297 for (
int col=0 ; col<nr ; col++) {
298 systeme.
set(ligne_courant, col+column_courant) = 0 ;
299 systeme.
set(ligne_courant, col+column_courant+nr) = (col+1)*(col+2)/double(2)*(1-2*(col%2)) ;
304 for (
int col=0 ; col<nr ; col++) {
305 systeme.
set(ligne_courant, col+column_courant) = 0 ;
306 systeme.
set(ligne_courant, col+column_courant+nr) = (col+1)*(col+2)/double(2)*(1-2*(col%2)) ;
307 systeme.
set(ligne_courant+1, col+column_courant) = 0 ;
308 systeme.
set(ligne_courant+1, col+column_courant+nr) = col*(col+1)*(col+2)*(col+3)/
double(12)*(2*(col%2)-1) ;
312 ligne_courant += nbr_cl ;
315 for (
int lig=0 ; lig<nr-1-nbr_cl ; lig++) {
316 for (
int col=0 ; col<nr ; col++) {
317 systeme.
set(lig+ligne_courant,col+column_courant) = - l_quant*(l_quant+1)*msn(lig,col) ;
318 sec_membre.
set(lig+ligne_courant) = 0 ;
319 systeme.
set(lig+ligne_courant,col+column_courant+nr)= mdn(lig,col) + 2*msn(lig,col) ;
325 ligne_courant += nr-1-nbr_cl ;
329 for (
int col=0 ; col<nr ; col++) {
332 systeme.
set(ligne_courant, col+column_courant) = 1 ;
333 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ;
336 systeme.
set(ligne_courant+1, col+column_courant) = 4*col*col/alpha ;
337 systeme.
set(ligne_courant+1, col+column_courant+nr) = 0 ;
340 systeme.
set(ligne_courant+1, col+column_courant) = (2*col+1)*(2*col+1)/alpha ;
341 systeme.
set(ligne_courant+1, col+column_courant+nr) = 0 ;
347 systeme.
set(ligne_courant, col+column_courant) = 1 ;
348 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ;
350 systeme.
set(ligne_courant+1, col+column_courant) = col*(col+3)/
double(2)/alpha ;
351 systeme.
set(ligne_courant+1, col+column_courant+nr) = 0 ;
357 for (
int col=0 ; col<nr ; col++) {
360 systeme.
set(ligne_courant, col+column_courant) = 0 ;
361 systeme.
set(ligne_courant, col+column_courant+nr) = 1 ;
364 systeme.
set(ligne_courant+1, col+column_courant) = 0 ;
365 systeme.
set(ligne_courant+1, col+column_courant+nr) = 4*col*col/alpha ;
368 systeme.
set(ligne_courant+1, col+column_courant) = 0 ;
369 systeme.
set(ligne_courant+1, col+column_courant+nr) = (2*col+1)*(2*col+1)/alpha ;
375 systeme.
set(ligne_courant, col+column_courant) = 0 ;
376 systeme.
set(ligne_courant, col+column_courant+nr) = 1 ;
378 systeme.
set(ligne_courant+1, col+column_courant) = 0 ;
379 systeme.
set(ligne_courant+1, col+column_courant+nr) = col*(col+3)/
double(2)/alpha ;
385 column_courant += 2*nr ;
390 for (
int l=1 ; l<nz-1 ; l++) {
393 alpha = (*mp_aff).get_alpha()[l] ;
394 beta = (*mp_aff).get_beta()[l] ;
401 for (
int col=0 ; col<nr ; col++) {
404 systeme.
set(ligne_courant, col+column_courant) = -1 ;
405 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ; }
407 systeme.
set(ligne_courant, col+column_courant) = 1 ;
408 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ; }
411 systeme.
set(ligne_courant+1, col+column_courant) = col*col/alpha ;
412 systeme.
set(ligne_courant+1, col+column_courant+nr) = 0 ; }
414 systeme.
set(ligne_courant+1, col+column_courant) = -col*col/alpha ;
415 systeme.
set(ligne_courant+1, col+column_courant+nr) = 0 ; }
420 for (
int col=0 ; col<nr ; col++) {
423 systeme.
set(ligne_courant, col+column_courant) = 0 ;
424 systeme.
set(ligne_courant, col+column_courant+nr) = -1 ; }
426 systeme.
set(ligne_courant, col+column_courant) = 0 ;
427 systeme.
set(ligne_courant, col+column_courant+nr) = 1 ; }
430 systeme.
set(ligne_courant+1, col+column_courant) = 0 ;
431 systeme.
set(ligne_courant+1, col+column_courant+nr) = col*col/alpha ; }
433 systeme.
set(ligne_courant+1, col+column_courant) = 0 ;
434 systeme.
set(ligne_courant+1, col+column_courant+nr) = -col*col/alpha ; }
444 for (
int i=0 ; i<nr ; i++)
446 Tbl xso(source_aux) ;
447 Tbl xxso(source_aux) ;
449 multx_1d(nr, &xxso.
t,
R_CHEB) ;
450 multx_1d(nr, &xxso.
t,
R_CHEB) ;
451 source_aux = beta*beta/alpha/alpha*source_aux+2*beta/alpha*xso+xxso ;
453 for (
int lig=0 ; lig<nr-2 ; lig++) {
454 for (
int col=0 ; col<nr ; col++) {
455 systeme.
set(lig+ligne_courant,col+column_courant) = mxds(lig,col) + mis(lig,col) ;
456 systeme.
set(lig+ligne_courant,col+column_courant+nr) = -mis(lig,col) ;
457 sec_membre.
set(lig+ligne_courant) = source_aux(lig) ;
458 systeme.
set(lig+ligne_courant+nr-2, col+column_courant) = -l_quant*(l_quant+1) ;
459 systeme.
set(lig+ligne_courant+nr-2, col+column_courant+nr) = mxds(lig,col) + 2*mis(lig,col) ;
460 sec_membre.
set(lig+ligne_courant+nr-2) = 0 ; }
464 ligne_courant += 2*nr-4 ;
466 for (
int col=0 ; col<nr ; col++) {
468 systeme.
set(ligne_courant, col+column_courant) = 1 ;
469 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ;
471 systeme.
set(ligne_courant+1, col+column_courant) = col*col/alpha ;
472 systeme.
set(ligne_courant+1, col+column_courant+nr) = 0 ;
474 systeme.
set(ligne_courant+2, col+column_courant) = 0 ;
475 systeme.
set(ligne_courant+2, col+column_courant+nr) = 1 ;
477 systeme.
set(ligne_courant+3, col+column_courant) = 0 ;
478 systeme.
set(ligne_courant+3, col+column_courant+nr) = col*col/alpha ;
481 column_courant += 2*nr ;
489 alpha = (*mp_aff).get_alpha()[nz-1] ;
490 beta = (*mp_aff).get_beta()[nz-1] ;
499 for (
int col=0 ; col<nr ; col++) {
502 systeme.
set(ligne_courant, col+column_courant) = -1 ;
503 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ; }
505 systeme.
set(ligne_courant, col+column_courant) = 1 ;
506 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ; }
509 systeme.
set(ligne_courant+1, col+column_courant) = -4*alpha*col*col ;
510 systeme.
set(ligne_courant+1, col+column_courant+nr) = 0 ; }
512 systeme.
set(ligne_courant+1, col+column_courant) = 4*alpha*col*col ;
513 systeme.
set(ligne_courant+1, col+column_courant+nr) = 0 ; }
516 systeme.
set(ligne_courant+2, col+column_courant) = -1 ;
517 systeme.
set(ligne_courant+2, col+column_courant+nr) = 0 ; }
519 systeme.
set(ligne_courant+2, col+column_courant) = 1 ;
520 systeme.
set(ligne_courant+2, col+column_courant+nr) = 0 ; }
523 systeme.
set(ligne_courant+3, col+column_courant) = -4*alpha*col*col ;
524 systeme.
set(ligne_courant+3, col+column_courant+nr) = 0 ; }
526 systeme.
set(ligne_courant+3, col+column_courant) = 4*alpha*col*col ;
527 systeme.
set(ligne_courant+3, col+column_courant+nr) = 0 ; }
538 for (
int col=0 ; col<nr ; col++) {
539 systeme.
set(ligne_courant, col+column_courant) = 1 ;
540 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ;
541 systeme.
set(ligne_courant+1, col+column_courant) = 0 ;
542 systeme.
set(ligne_courant+1, col+column_courant+nr) = 1 ; }
547 for (
int col=0 ; col<nr ; col++) {
548 systeme.
set(ligne_courant, col+column_courant) = 1 ;
549 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ;
550 systeme.
set(ligne_courant+1, col+column_courant) = 0 ;
551 systeme.
set(ligne_courant+1, col+column_courant+nr) = 1; }
553 for (
int col=0 ; col<nr ; col++) {
554 systeme.
set(ligne_courant+2, col+column_courant) = -4*alpha*col*col ;
555 systeme.
set(ligne_courant+2, col+column_courant+nr) = 0 ;
556 systeme.
set(ligne_courant+3, col+column_courant) = 0 ;
557 systeme.
set(ligne_courant+3, col+column_courant+nr) = -4*alpha*col*col ; }
564 for (
int col=0 ; col<nr ; col++) {
565 systeme.
set(ligne_courant, col+column_courant) = 1 ;
566 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ;
567 systeme.
set(ligne_courant+1, col+column_courant) = 0 ;
568 systeme.
set(ligne_courant+1, col+column_courant+nr) = 1 ; }
575 for (
int col=0 ; col<nr ; col++) {
576 systeme.
set(ligne_courant, col+column_courant) = 1 ;
577 systeme.
set(ligne_courant, col+column_courant+nr) = 0 ;
578 systeme.
set(ligne_courant+1, col+column_courant) = 0 ;
579 systeme.
set(ligne_courant+1, col+column_courant+1) = 1 ; }
583 cout <<
"Unknown dzpuis in vector_divfree_A ..." << endl ;
587 ligne_courant += nbr_cl ;
593 indic = alpha*alpha ;
603 for (
int lig=0 ; lig<nr-1-nbr_cl ; lig++) {
604 for (
int col=0 ; col<nr ; col++) {
605 systeme.
set(lig+ligne_courant,col+column_courant) = mdzec(lig,col)+mszec(lig,col) ;
606 systeme.
set(lig+ligne_courant,col+column_courant+nr) = -mszec(lig,col) ;
608 systeme.
set(lig+ligne_courant+nr-1-nbr_cl,col+column_courant)= - l_quant*(l_quant+1)*mszec(lig,col) ;
609 systeme.
set(lig+ligne_courant+nr-1-nbr_cl,col+column_courant+nr) = mdzec(lig,col)+2*mszec(lig,col) ;
610 sec_membre.
set(lig+ligne_courant+nr-1-nbr_cl) = 0 ; }
621 for (
int l=0 ; l<nz ; l++) {
623 for (
int i=0 ; i<nr ; i++) {
624 meta.
set(l,k,j,i) = soluce(conte) ;
625 mvr.
set(l,k,i,j) = soluce(conte+nr) ;