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) ;
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
void ylm_i()
Inverse of ylm()
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
void ylm()
Computes the coefficients of *this.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
double & set(int i)
Read/write of a particular element (index i) (1D case)
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Tensor field of valence 0 (or component of a tensorial field).
Class for the elementary differential operator (see the base class Diff ).
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
void annule_hard()
Sets the Scalar to zero in a hard way.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Class for the elementary differential operator Identity (see the base class Diff ).
int get_n_int() const
Returns the number of stored int 's addresses.
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
int get_dzpuis() const
Returns dzpuis.
double * t
The array of double.
Mtbl * c
Values of the function at the points of the multi-grid.
int get_nzone() const
Returns the number of domains.
void annule_hard()
Sets the logical state to ETATQCQ (undefined state).
double & set(int j, int i)
Read/write of a particuliar element.
void set_band(int up, int low) const
Calculate the band storage of *std.
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Bases of the spectral expansions.
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Class for the elementary differential operator (see the base class Diff ).
void mult_x()
The basis is transformed as with a multiplication by .
Coefficients storage for the multi-domain spectral method.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Valeur & set_spectral_va()
Returns va (read/write version)
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 )...
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
const Map & get_mp() const
Returns the mapping.
void annule_hard()
Sets the Tbl to zero in a hard way.
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
const Valeur & get_spectral_va() const
Returns va (read only version)
void sol_Dirac_A_tau(const Scalar &aaa, Scalar &eta, Scalar &vr, const Param *par_bc=0x0) const
Solves via a tau method a system of two-coupled first-order PDEs obtained from the divergence-free co...
#define R_CHEB
base de Chebychev ordinaire (fin)