27 char scalar_match_tau_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_match_tau.C,v 1.6 2014/10/13 08:53:46 j_novak Exp $" ;
69 assert(mp_aff != 0x0) ;
75 if (
etat == ETATZERO) {
76 if (mu.get_etat() == ETATZERO)
84 bool ced = (mgrid.
get_type_r(nzm1) == UNSURR) ;
86 int domain_bc = par_bc.
get_int(0) ;
87 bool bc_ced = ((ced) && (domain_bc == nzm1)) ;
89 int n_conditions = par_bc.
get_int(1) ;
90 assert ((n_conditions==1)||(n_conditions==2)) ;
91 bool derivative = (n_conditions == 2) ;
92 int dl = 0 ;
int l_min = 0 ;
int excised_i =0;
100 bool isexcised = (excised_i==1);
102 int nt = mgrid.
get_nt(0) ;
103 int np = mgrid.
get_np(0) ;
112 if (mu.get_etat() == ETATZERO)
114 int system01_size =0;
116 if (isexcised ==
false){
124 for (
int lz=1; lz<=domain_bc; lz++) {
125 system01_size += n_conditions ;
126 system_size += n_conditions ;
128 assert (
etat != ETATNONDEF) ;
130 bool need_matrices = true ;
133 need_matrices =
false ;
135 Matrice system_l0(system01_size, system01_size) ;
136 Matrice system_l1(system01_size, system01_size) ;
137 Matrice system_even(system_size, system_size) ;
138 Matrice system_odd(system_size, system_size) ;
145 int index = 0 ;
int index01 = 0 ;
148 if (isexcised ==
false){
149 system_even.
set(index, index) = 1. ;
150 system_even.
set(index, index + 1) = -1. ;
151 system_odd.
set(index, index) = -(2.*nr - 5.)/alpha ;
152 system_odd.
set(index, index+1) = (2.*nr - 3.)/alpha ;
154 if (domain_bc == 0) {
155 system_l0.
set(index01, index01) = c_val + c_der*4.*(nr-1)*(nr-1)/alpha ;
156 system_l1.
set(index01, index01) = c_val + c_der*(2*nr-3)*(2*nr-3)/alpha ;
158 system_even.
set(index, index-1) = c_val + c_der*4.*(nr-2)*(nr-2)/alpha ;
159 system_even.
set(index, index) = c_val + c_der*4.*(nr-1)*(nr-1)/alpha ;
160 system_odd.
set(index, index-1) = c_val + c_der*(2*nr-5)*(2*nr-5)/alpha ;
161 system_odd.
set(index, index) = c_val + c_der*(2*nr-3)*(2*nr-3)/alpha ;
165 system_l0.
set(index01, index01) = 1. ;
166 system_l1.
set(index01, index01) = 1. ;
167 system_even.
set(index, index-1) = 1. ;
168 system_even.
set(index, index) = 1. ;
169 system_odd.
set(index, index-1) = 1. ;
170 system_odd.
set(index, index) = 1. ;
172 system_l0.
set(index01+1, index01) = 4*(nr-1)*(nr-1)/alpha ;
173 system_l1.
set(index01+1, index01) = (2*nr-3)*(2*nr-3)/alpha ;
174 system_even.
set(index+1, index-1) = 4*(nr-2)*(nr-2)/alpha ;
175 system_even.
set(index+1, index) = 4*(nr-1)*(nr-1)/alpha ;
176 system_odd.
set(index+1, index-1) = (2*nr-5)*(2*nr-5)/alpha ;
177 system_odd.
set(index+1, index) = (2*nr-3)*(2*nr-3)/alpha ;
184 if ((1 == domain_bc)&&(bc_ced))
185 alpha = -0.25/alpha ;
186 system_l0.
set(index01, index01+1) = 1. ;
187 system_l0.
set(index01, index01+2) = -1. ;
188 system_l1.
set(index01, index01+1) = 1. ;
189 system_l1.
set(index01, index01+2) = -1. ;
191 system_even.
set(index, index+1) = 1. ;
192 system_even.
set(index, index+2) = -1. ;
193 system_odd.
set(index, index+1) = 1. ;
194 system_odd.
set(index, index+2) = -1. ;
197 system_l0.
set(index01, index01) = -(nr-2)*(nr-2)/alpha ;
198 system_l0.
set(index01, index01+1) = (nr-1)*(nr-1)/alpha ;
199 system_l1.
set(index01, index01) = -(nr-2)*(nr-2)/alpha ;
200 system_l1.
set(index01, index01+1) = (nr-1)*(nr-1)/alpha ;
202 system_even.
set(index, index) = -(nr-2)*(nr-2)/alpha ;
203 system_even.
set(index, index+1) = (nr-1)*(nr-1)/alpha ;
204 system_odd.
set(index, index) = -(nr-2)*(nr-2)/alpha ;
205 system_odd.
set(index, index+1) = (nr-1)*(nr-1)/alpha ;
209 if (1 == domain_bc) {
210 system_l0.
set(index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
211 system_l0.
set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
212 system_l1.
set(index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
213 system_l1.
set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
215 system_even.
set(index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
216 system_even.
set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
217 system_odd.
set(index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
218 system_odd.
set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
222 system_l0.
set(index01, index01-1) = 1. ;
223 system_l0.
set(index01, index01) = 1. ;
224 system_l1.
set(index01, index01-1) = 1. ;
225 system_l1.
set(index01, index01) = 1. ;
226 system_even.
set(index, index-1) = 1. ;
227 system_even.
set(index, index) = 1. ;
228 system_odd.
set(index, index-1) = 1. ;
229 system_odd.
set(index, index) = 1. ;
231 system_l0.
set(index01+1, index01-1) = (nr-2)*(nr-2)/alpha ;
232 system_l0.
set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
233 system_l1.
set(index01+1, index01-1) = (nr-2)*(nr-2)/alpha ;
234 system_l1.
set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
235 system_even.
set(index+1, index-1) = (nr-2)*(nr-2)/alpha ;
236 system_even.
set(index+1, index) = (nr-1)*(nr-1)/alpha ;
237 system_odd.
set(index+1, index-1) = (nr-2)*(nr-2)/alpha ;
238 system_odd.
set(index+1, index) = (nr-1)*(nr-1)/alpha ;
246 if ((1 == domain_bc)&&(bc_ced))
247 alpha = -0.25/alpha ;
248 if (1 == domain_bc) {
249 system_l0.
set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
250 system_l1.
set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
252 system_even.
set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
253 system_odd.
set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
257 system_l0.
set(index01, index01) = 1. ;
258 system_l1.
set(index01, index01) = 1. ;
259 system_even.
set(index, index) = 1. ;
260 system_odd.
set(index, index) = 1. ;
262 system_l0.
set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
263 system_l1.
set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
264 system_even.
set(index+1, index) = (nr-1)*(nr-1)/alpha ;
265 system_odd.
set(index+1, index) = (nr-1)*(nr-1)/alpha ;
270 for (
int lz=2; lz<=domain_bc; lz++) {
273 if ((lz == domain_bc)&&(bc_ced))
274 alpha = -0.25/alpha ;
275 system_l0.
set(index01, index01+1) = 1. ;
276 system_l0.
set(index01, index01+2) = -1. ;
277 system_l1.
set(index01, index01+1) = 1. ;
278 system_l1.
set(index01, index01+2) = -1. ;
280 system_even.
set(index, index+1) = 1. ;
281 system_even.
set(index, index+2) = -1. ;
282 system_odd.
set(index, index+1) = 1. ;
283 system_odd.
set(index, index+2) = -1. ;
286 system_l0.
set(index01, index01) = -(nr-2)*(nr-2)/alpha ;
287 system_l0.
set(index01, index01+1) = (nr-1)*(nr-1)/alpha ;
288 system_l1.
set(index01, index01) = -(nr-2)*(nr-2)/alpha ;
289 system_l1.
set(index01, index01+1) = (nr-1)*(nr-1)/alpha ;
291 system_even.
set(index, index) = -(nr-2)*(nr-2)/alpha ;
292 system_even.
set(index, index+1) = (nr-1)*(nr-1)/alpha ;
293 system_odd.
set(index, index) = -(nr-2)*(nr-2)/alpha ;
294 system_odd.
set(index, index+1) = (nr-1)*(nr-1)/alpha ;
298 if (lz == domain_bc) {
299 system_l0.
set(index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
300 system_l0.
set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
301 system_l1.
set(index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
302 system_l1.
set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
304 system_even.
set(index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
305 system_even.
set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
306 system_odd.
set(index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
307 system_odd.
set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
311 system_l0.
set(index01, index01-1) = 1. ;
312 system_l0.
set(index01, index01) = 1. ;
313 system_l1.
set(index01, index01-1) = 1. ;
314 system_l1.
set(index01, index01) = 1. ;
315 system_even.
set(index, index-1) = 1. ;
316 system_even.
set(index, index) = 1. ;
317 system_odd.
set(index, index-1) = 1. ;
318 system_odd.
set(index, index) = 1. ;
320 system_l0.
set(index01+1, index01-1) = (nr-2)*(nr-2)/alpha ;
321 system_l0.
set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
322 system_l1.
set(index01+1, index01-1) = (nr-2)*(nr-2)/alpha ;
323 system_l1.
set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
324 system_even.
set(index+1, index-1) = (nr-2)*(nr-2)/alpha ;
325 system_even.
set(index+1, index) = (nr-1)*(nr-1)/alpha ;
326 system_odd.
set(index+1, index-1) = (nr-2)*(nr-2)/alpha ;
327 system_odd.
set(index+1, index) = (nr-1)*(nr-1)/alpha ;
332 assert(index01 == system01_size) ;
333 assert(index == system_size) ;
338 if (par_mat != 0x0) {
352 const Matrice& sys_even = (need_matrices ? system_even :
354 const Matrice& sys_odd = (need_matrices ? system_odd :
365 int m_q, l_q, base_r ;
366 for (
int k=0; k<np+2; k++) {
367 for (
int j=0; j<nt; j++) {
369 if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q >= l_min)) {
371 int sys_size = (l_q < 2 ? system01_size : system_size) ;
372 int nl = (l_q < 2 ? 1 : 2) ;
377 int nr = mgrid.
get_nr(0) ;
380 if (isexcised==
false){
383 double val_c = 0.; pari = 1 ;
384 for (
int i=0; i<nr-nl; i++) {
385 val_c += pari*coef(0, k, j, i) ;
388 rhs.
set(index) = val_c ;
392 double der_c = 0.; pari = 1 ;
393 for (
int i=0; i<nr-nl-1; i++) {
394 der_c += pari*(2*i+1)*coef(0, k, j, i) ;
397 rhs.
set(index) = der_c / alpha ;
402 for (
int i=0; i<nr-nl; i++) {
403 val_b += coef(0, k, j, i) ;
404 der_b += 4*i*i*coef(0, k, j, i) ;
409 for (
int i=0; i<nr-nl-1; i++) {
410 val_b += coef(0, k, j, i) ;
411 der_b += (2*i+1)*(2*i+1)*coef(0, k, j, i) ;
415 rhs.
set(index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha ;
419 rhs.
set(index) = -val_b ;
421 rhs.
set(index+1) = -der_b/alpha ;
425 if ((1 == domain_bc)&&(bc_ced))
426 alpha = -0.25/alpha ;
428 int nr_lim = nr - n_conditions ;
429 val_b = 0 ; pari = 1 ;
430 for (
int i=0; i<nr_lim; i++) {
431 val_b += pari*coef(1, k, j, i) ;
434 rhs.
set(index) += val_b ;
437 der_b = 0 ; pari = -1 ;
438 for (
int i=0; i<nr_lim; i++) {
439 der_b += pari*i*i*coef(1, k, j, i) ;
442 rhs.
set(index) += der_b/alpha ;
446 for (
int i=0; i<nr_lim; i++)
447 val_b += coef(1, k, j, i) ;
449 for (
int i=0; i<nr_lim; i++) {
450 der_b += i*i*coef(1, k, j, i) ;
453 rhs.
set(index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha ;
457 rhs.
set(index) = -val_b ;
458 if (derivative) rhs.
set(index+1) = -der_b / alpha ;
464 if ((1 == domain_bc)&&(bc_ced))
465 alpha = -0.25/alpha ;
467 int nr_lim = nr - 1 ;
469 for (
int i=0; i<nr_lim; i++)
470 val_b += coef(1, k, j, i) ;
472 for (
int i=0; i<nr_lim; i++) {
473 der_b += i*i*coef(1, k, j, i) ;
476 rhs.
set(index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha ;
480 rhs.
set(index) = -val_b ;
481 if (derivative) rhs.
set(index+1) = -der_b / alpha ;
487 for (
int lz=2; lz<=domain_bc; lz++) {
489 if ((lz == domain_bc)&&(bc_ced))
490 alpha = -0.25/alpha ;
492 int nr_lim = nr - n_conditions ;
493 val_b = 0 ; pari = 1 ;
494 for (
int i=0; i<nr_lim; i++) {
495 val_b += pari*coef(lz, k, j, i) ;
498 rhs.
set(index) += val_b ;
501 der_b = 0 ; pari = -1 ;
502 for (
int i=0; i<nr_lim; i++) {
503 der_b += pari*i*i*coef(lz, k, j, i) ;
506 rhs.
set(index) += der_b/alpha ;
510 for (
int i=0; i<nr_lim; i++)
511 val_b += coef(lz, k, j, i) ;
513 for (
int i=0; i<nr_lim; i++) {
514 der_b += i*i*coef(lz, k, j, i) ;
517 rhs.
set(index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha ;
521 rhs.
set(index) = -val_b ;
522 if (derivative) rhs.
set(index+1) = -der_b / alpha ;
535 solut = (l_q%2 == 0 ? sys_even.
inverse(rhs) :
539 int dec = (base_r ==
R_CHEBP ? 0 : 1) ;
541 if (isexcised==
false){
543 coef.
set(0, k, j, nr-2-dec) = solut(index) ;
546 coef.
set(0, k, j, nr-1-dec) = solut(index) ;
549 coef.
set(0, k, j, nr-1) = 0 ;
553 for (nl=1; nl<=n_conditions; nl++) {
554 int ii = n_conditions - nl + 1 ;
555 coef.
set(1, k, j, nr-ii) = solut(index) ;
561 coef.
set(1,k,j,nr-1)=solut(index);
564 for (
int lz=2; lz<=domain_bc; lz++) {
566 for (nl=1; nl<=n_conditions; nl++) {
567 int ii = n_conditions - nl + 1 ;
568 coef.
set(lz, k, j, nr-ii) = solut(index) ;
574 for (
int lz=0; lz<=domain_bc; lz++)
575 for (
int i=0; i<mgrid.
get_nr(lz); i++)
576 coef.
set(lz, k, j, i) = 0 ;
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
double & get_double_mod(int position=0) const
Returns the reference of a stored modifiable double .
const double * get_alpha() const
Returns the pointer on the array alpha.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
void coef() const
Computes the coeffcients of *this.
int get_n_matrice_mod() const
Returns the number of modifiable Matrice 's addresses in the list.
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)
void add_matrice_mod(Matrice &ti, int position=0)
Adds the address of a new modifiable Matrice to the list.
void match_tau(Param &par_bc, Param *par_mat=0x0)
Method for matching accross domains and imposing boundary condition.
void annule_hard()
Sets the Scalar to zero in a hard way.
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.
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define R_CHEBP
base de Cheb. paire (rare) seulement
Mtbl * c
Values of the function at the points of the multi-grid.
Tbl & get_tbl_mod(int position=0) const
Returns the reference of a modifiable Tbl stored in the list.
int get_n_double_mod() const
Returns the number of stored modifiable double 's addresses.
int get_nzone() const
Returns the number of domains.
Valeur va
The numerical value of the Scalar.
void annule_hard()
Sets the logical state to ETATQCQ (undefined state).
double & set(int j, int i)
Read/write of a particuliar element.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
int get_n_tbl_mod() const
Returns the number of modifiable Tbl 's addresses in the list.
Bases of the spectral expansions.
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary)...
Coefficients storage for the multi-domain spectral method.
Matrice & get_matrice_mod(int position=0) const
Returns the reference of a modifiable Matrice stored in the list.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
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.
void annule_hard()
Sets the Tbl to zero in a hard way.