LORENE
scalar_match_tau.C
1 /*
2  * Function to perform a matching across domains + imposition of boundary conditions
3  * at the outer domain and regularity condition at the center.
4  */
5 
6 /*
7  * Copyright (c) 2008 Jerome Novak
8  * 2011 Nicolas Vasset
9  *
10  * This file is part of LORENE.
11  *
12  * LORENE is free software; you can redistribute it and/or modify
13  * it under the terms of the GNU General Public License version 2
14  * as published by the Free Software Foundation.
15  *
16  * LORENE is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  * GNU General Public License for more details.
20  *
21  * You should have received a copy of the GNU General Public License
22  * along with LORENE; if not, write to the Free Software
23  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24  *
25  */
26 
27 
28 
29 /*
30  * $Id: scalar_match_tau.C,v 1.9 2024/06/25 14:03:31 j_novak Exp $
31  * $Log: scalar_match_tau.C,v $
32  * Revision 1.9 2024/06/25 14:03:31 j_novak
33  * Minor changes
34  *
35  * Revision 1.8 2024/01/26 17:44:26 g_servignat
36  * Updated the Pseudopolytrope_1D class to be consistent with the paper (i.e. with a GPP in the middle)
37  *
38  * Revision 1.7 2016/12/05 16:18:18 j_novak
39  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
40  *
41  * Revision 1.6 2014/10/13 08:53:46 j_novak
42  * Lorene classes and functions now belong to the namespace Lorene.
43  *
44  * Revision 1.5 2014/10/06 15:16:16 j_novak
45  * Modified #include directives to use c++ syntax.
46  *
47  * Revision 1.4 2012/05/11 14:11:24 j_novak
48  * Modifications to deal with the accretion of a scalar field
49  *
50  * Revision 1.3 2012/02/06 12:59:07 j_novak
51  * Correction of some errors.
52  *
53  * Revision 1.2 2008/08/20 13:23:43 j_novak
54  * The shift in quantum number l (e.g. for \tilde{B}) is now taken into account.
55  *
56  * Revision 1.1 2008/05/24 15:05:22 j_novak
57  * New method Scalar::match_tau to match the output of an explicit time-marching scheme with the tau method.
58  *
59  *
60  * $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_match_tau.C,v 1.9 2024/06/25 14:03:31 j_novak Exp $
61  *
62  */
63 
64 // C headers
65 #include <cassert>
66 #include <cmath>
67 
68 // Lorene headers
69 #include "tensor.h"
70 #include "matrice.h"
71 #include "proto.h"
72 #include "param.h"
73 
74 namespace Lorene {
75 void Scalar::match_tau(Param& par_bc, Param* par_mat) {
76 
77  const Map_af* mp_aff = dynamic_cast<const Map_af*>(mp) ;
78  assert(mp_aff != 0x0) ; //Only affine mapping for the moment
79 
80  const Mg3d& mgrid = *mp_aff->get_mg() ;
81  assert(mgrid.get_type_r(0) == RARE) ;
82  assert (par_bc.get_n_tbl_mod() > 0) ;
83  Tbl mu = par_bc.get_tbl_mod(0) ;
84  if (etat == ETATZERO) {
85  if (mu.get_etat() == ETATZERO)
86  return ;
87  else
88  annule_hard() ;
89  }
90 
91  int nz = mgrid.get_nzone() ;
92  int nzm1 = nz - 1 ;
93  bool ced = (mgrid.get_type_r(nzm1) == UNSURR) ;
94  assert(par_bc.get_n_int() >= 2);
95  int domain_bc = par_bc.get_int(0) ;
96  bool bc_ced = ((ced) && (domain_bc == nzm1)) ;
97 
98  int n_conditions = par_bc.get_int(1) ;
99  assert ((n_conditions==1)||(n_conditions==2)) ;
100  bool derivative = (n_conditions == 2) ;
101  int dl = 0 ; int l_min = 0 ; int excised_i =0;
102  if (par_bc.get_n_int() > 2) {
103  dl = par_bc.get_int(2) ;
104  l_min = par_bc.get_int(3) ;
105  if(par_bc.get_n_int() >4){
106  excised_i = par_bc.get_int(4);
107  }
108  }
109  bool isexcised = (excised_i==1);
110 
111  int nt = mgrid.get_nt(0) ;
112  int np = mgrid.get_np(0) ;
113  assert (par_bc.get_n_double_mod() == 2) ;
114  double c_val = par_bc.get_double_mod(0) ;
115  double c_der = par_bc.get_double_mod(1) ;
116  if (bc_ced) {
117  c_val = 1 ;
118  c_der = 0 ;
119  mu.annule_hard() ;
120  }
121  if (mu.get_etat() == ETATZERO)
122  mu.annule_hard() ;
123  int system01_size =0;
124  int system_size =0;
125  if (isexcised ==false){
126  system01_size = 1 ;
127  system_size = 2 ;
128  }
129  else{
130  system01_size = -1;
131  system_size = -1;
132  }
133  for (int lz=1; lz<=domain_bc; lz++) {
134  system01_size += n_conditions ;
135  system_size += n_conditions ;
136  }
137  assert (etat != ETATNONDEF) ;
138 
139  bool need_matrices = true ;
140  if (par_mat != 0x0)
141  if (par_mat->get_n_matrice_mod() == 4)
142  need_matrices = false ;
143 
144  Matrice system_l0(system01_size, system01_size) ;
145  Matrice system_l1(system01_size, system01_size) ;
146  Matrice system_even(system_size, system_size) ;
147  Matrice system_odd(system_size, system_size) ;
148 
149  if (need_matrices) {
150  system_l0.annule_hard() ;
151  system_l1.annule_hard() ;
152  system_even.annule_hard() ;
153  system_odd.annule_hard() ;
154  int index = 0 ; int index01 = 0 ;
155  int nr = mgrid.get_nr(0);
156  double alpha = mp_aff->get_alpha()[0] ;
157  if (isexcised == false){
158  system_even.set(index, index) = 1. ;
159  system_even.set(index, index + 1) = -1. ; //C_{N-1} - C_N = \sum (-1)^i C_i
160  system_odd.set(index, index) = -(2.*nr - 5.)/alpha ;
161  system_odd.set(index, index+1) = (2.*nr - 3.)/alpha ;
162  index++ ;
163  if (domain_bc == 0) {
164  system_l0.set(index01, index01) = c_val + c_der*4.*(nr-1)*(nr-1)/alpha ;
165  system_l1.set(index01, index01) = c_val + c_der*(2*nr-3)*(2*nr-3)/alpha ;
166  index01++ ;
167  system_even.set(index, index-1) = c_val + c_der*4.*(nr-2)*(nr-2)/alpha ;
168  system_even.set(index, index) = c_val + c_der*4.*(nr-1)*(nr-1)/alpha ;
169  system_odd.set(index, index-1) = c_val + c_der*(2*nr-5)*(2*nr-5)/alpha ;
170  system_odd.set(index, index) = c_val + c_der*(2*nr-3)*(2*nr-3)/alpha ;
171  index++ ;
172  }
173  else {
174  system_l0.set(index01, index01) = 1. ;
175  system_l1.set(index01, index01) = 1. ;
176  system_even.set(index, index-1) = 1. ;
177  system_even.set(index, index) = 1. ;
178  system_odd.set(index, index-1) = 1. ;
179  system_odd.set(index, index) = 1. ;
180  if (derivative) {
181  system_l0.set(index01+1, index01) = 4*(nr-1)*(nr-1)/alpha ;
182  system_l1.set(index01+1, index01) = (2*nr-3)*(2*nr-3)/alpha ;
183  system_even.set(index+1, index-1) = 4*(nr-2)*(nr-2)/alpha ;
184  system_even.set(index+1, index) = 4*(nr-1)*(nr-1)/alpha ;
185  system_odd.set(index+1, index-1) = (2*nr-5)*(2*nr-5)/alpha ;
186  system_odd.set(index+1, index) = (2*nr-3)*(2*nr-3)/alpha ;
187  }
188  }
189  if (domain_bc > 0) {
190  // Do things for lz=1;
191  nr = mgrid.get_nr(1) ;
192  alpha = mp_aff->get_alpha()[1] ;
193  if ((1 == domain_bc)&&(bc_ced))
194  alpha = -0.25/alpha ;
195  system_l0.set(index01, index01+1) = 1. ;
196  system_l0.set(index01, index01+2) = -1. ;
197  system_l1.set(index01, index01+1) = 1. ;
198  system_l1.set(index01, index01+2) = -1. ;
199  index01++ ;
200  system_even.set(index, index+1) = 1. ;
201  system_even.set(index, index+2) = -1. ;
202  system_odd.set(index, index+1) = 1. ;
203  system_odd.set(index, index+2) = -1. ;
204  index++ ;
205  if (derivative) {
206  system_l0.set(index01, index01) = -(nr-2)*(nr-2)/alpha ;
207  system_l0.set(index01, index01+1) = (nr-1)*(nr-1)/alpha ;
208  system_l1.set(index01, index01) = -(nr-2)*(nr-2)/alpha ;
209  system_l1.set(index01, index01+1) = (nr-1)*(nr-1)/alpha ;
210  index01++ ;
211  system_even.set(index, index) = -(nr-2)*(nr-2)/alpha ;
212  system_even.set(index, index+1) = (nr-1)*(nr-1)/alpha ;
213  system_odd.set(index, index) = -(nr-2)*(nr-2)/alpha ;
214  system_odd.set(index, index+1) = (nr-1)*(nr-1)/alpha ;
215  index++ ;
216  }
217 
218  if (1 == domain_bc) {
219  system_l0.set(index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
220  system_l0.set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
221  system_l1.set(index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
222  system_l1.set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
223  index01++ ;
224  system_even.set(index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
225  system_even.set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
226  system_odd.set(index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
227  system_odd.set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
228  index++ ;
229  }
230  else {
231  system_l0.set(index01, index01-1) = 1. ;
232  system_l0.set(index01, index01) = 1. ;
233  system_l1.set(index01, index01-1) = 1. ;
234  system_l1.set(index01, index01) = 1. ;
235  system_even.set(index, index-1) = 1. ;
236  system_even.set(index, index) = 1. ;
237  system_odd.set(index, index-1) = 1. ;
238  system_odd.set(index, index) = 1. ;
239  if (derivative) {
240  system_l0.set(index01+1, index01-1) = (nr-2)*(nr-2)/alpha ;
241  system_l0.set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
242  system_l1.set(index01+1, index01-1) = (nr-2)*(nr-2)/alpha ;
243  system_l1.set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
244  system_even.set(index+1, index-1) = (nr-2)*(nr-2)/alpha ;
245  system_even.set(index+1, index) = (nr-1)*(nr-1)/alpha ;
246  system_odd.set(index+1, index-1) = (nr-2)*(nr-2)/alpha ;
247  system_odd.set(index+1, index) = (nr-1)*(nr-1)/alpha ;
248  }
249  }
250  }
251  }
252  else {
253  nr = mgrid.get_nr(1) ;
254  alpha = mp_aff->get_alpha()[1] ;
255  if ((1 == domain_bc)&&(bc_ced))
256  alpha = -0.25/alpha ;
257  if (1 == domain_bc) {
258  system_l0.set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
259  system_l1.set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
260  index01++ ;
261  system_even.set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
262  system_odd.set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
263  index++ ;
264  }
265  else {
266  system_l0.set(index01, index01) = 1. ;
267  system_l1.set(index01, index01) = 1. ;
268  system_even.set(index, index) = 1. ;
269  system_odd.set(index, index) = 1. ;
270  if (derivative) {
271  system_l0.set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
272  system_l1.set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
273  system_even.set(index+1, index) = (nr-1)*(nr-1)/alpha ;
274  system_odd.set(index+1, index) = (nr-1)*(nr-1)/alpha ;
275  }
276 
277  }
278  }
279  for (int lz=2; lz<=domain_bc; lz++) {
280  nr = mgrid.get_nr(lz) ;
281  alpha = mp_aff->get_alpha()[lz] ;
282  if ((lz == domain_bc)&&(bc_ced))
283  alpha = -0.25/alpha ;
284  system_l0.set(index01, index01+1) = 1. ;
285  system_l0.set(index01, index01+2) = -1. ;
286  system_l1.set(index01, index01+1) = 1. ;
287  system_l1.set(index01, index01+2) = -1. ;
288  index01++ ;
289  system_even.set(index, index+1) = 1. ;
290  system_even.set(index, index+2) = -1. ;
291  system_odd.set(index, index+1) = 1. ;
292  system_odd.set(index, index+2) = -1. ;
293  index++ ;
294  if (derivative) {
295  system_l0.set(index01, index01) = -(nr-2)*(nr-2)/alpha ;
296  system_l0.set(index01, index01+1) = (nr-1)*(nr-1)/alpha ;
297  system_l1.set(index01, index01) = -(nr-2)*(nr-2)/alpha ;
298  system_l1.set(index01, index01+1) = (nr-1)*(nr-1)/alpha ;
299  index01++ ;
300  system_even.set(index, index) = -(nr-2)*(nr-2)/alpha ;
301  system_even.set(index, index+1) = (nr-1)*(nr-1)/alpha ;
302  system_odd.set(index, index) = -(nr-2)*(nr-2)/alpha ;
303  system_odd.set(index, index+1) = (nr-1)*(nr-1)/alpha ;
304  index++ ;
305  }
306 
307  if (lz == domain_bc) {
308  system_l0.set(index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
309  system_l0.set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
310  system_l1.set(index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
311  system_l1.set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
312  index01++ ;
313  system_even.set(index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
314  system_even.set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
315  system_odd.set(index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
316  system_odd.set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
317  index++ ;
318  }
319  else {
320  system_l0.set(index01, index01-1) = 1. ;
321  system_l0.set(index01, index01) = 1. ;
322  system_l1.set(index01, index01-1) = 1. ;
323  system_l1.set(index01, index01) = 1. ;
324  system_even.set(index, index-1) = 1. ;
325  system_even.set(index, index) = 1. ;
326  system_odd.set(index, index-1) = 1. ;
327  system_odd.set(index, index) = 1. ;
328  if (derivative) {
329  system_l0.set(index01+1, index01-1) = (nr-2)*(nr-2)/alpha ;
330  system_l0.set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
331  system_l1.set(index01+1, index01-1) = (nr-2)*(nr-2)/alpha ;
332  system_l1.set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
333  system_even.set(index+1, index-1) = (nr-2)*(nr-2)/alpha ;
334  system_even.set(index+1, index) = (nr-1)*(nr-1)/alpha ;
335  system_odd.set(index+1, index-1) = (nr-2)*(nr-2)/alpha ;
336  system_odd.set(index+1, index) = (nr-1)*(nr-1)/alpha ;
337  }
338  }
339  } // End of loop on lz
340 
341  assert(index01 == system01_size) ;
342  assert(index == system_size) ;
343  system_l0.set_lu() ;
344  system_l1.set_lu() ;
345  system_even.set_lu() ;
346  system_odd.set_lu() ;
347  if (par_mat != 0x0) {
348  Matrice* pmat0 = new Matrice(system_l0) ;
349  Matrice* pmat1 = new Matrice(system_l1) ;
350  Matrice* pmat2 = new Matrice(system_even) ;
351  Matrice* pmat3 = new Matrice(system_odd) ;
352  par_mat->add_matrice_mod(*pmat0, 0) ;
353  par_mat->add_matrice_mod(*pmat1, 1) ;
354  par_mat->add_matrice_mod(*pmat2, 2) ;
355  par_mat->add_matrice_mod(*pmat3, 3) ;
356  }
357  }// End of if (need_matrices) ...
358 
359  const Matrice& sys_l0 = (need_matrices ? system_l0 : par_mat->get_matrice_mod(0) ) ;
360  const Matrice& sys_l1 = (need_matrices ? system_l1 : par_mat->get_matrice_mod(1) ) ;
361  const Matrice& sys_even = (need_matrices ? system_even :
362  par_mat->get_matrice_mod(2) ) ;
363  const Matrice& sys_odd = (need_matrices ? system_odd :
364  par_mat->get_matrice_mod(3) ) ;
365 
366  va.coef() ;
367  va.ylm() ;
368  const Base_val& base = get_spectral_base() ;
369  Mtbl_cf& coef = *va.c_cf ;
370  if (va.c != 0x0) {
371  delete va.c ;
372  va.c = 0x0 ;
373  }
374  int m_q, l_q, base_r ;
375  for (int k=0; k<np+2; k++) {
376  for (int j=0; j<nt; j++) {
377  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;//#0 here as domain index
378  if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q >= l_min)) {
379  l_q += dl ;
380  int sys_size = (l_q < 2 ? system01_size : system_size) ;
381  int nl = (l_q < 2 ? 1 : 2) ;
382  Tbl rhs(sys_size) ; rhs.annule_hard() ;
383  int index = 0 ;
384  int pari = 1 ;
385  double alpha = mp_aff->get_alpha()[0] ;
386  int nr = mgrid.get_nr(0) ;
387  double val_b = 0. ;
388  double der_b =0. ;
389  if (isexcised==false){
390  if (l_q>=2) {
391  if (base_r == R_CHEBP) {
392  double val_c = 0.; pari = 1 ;
393  for (int i=0; i<nr-nl; i++) {
394  val_c += pari*coef(0, k, j, i) ;
395  pari = -pari ;
396  }
397  rhs.set(index) = val_c ;
398  }
399  else {
400  assert( base_r == R_CHEBI) ;
401  double der_c = 0.; pari = 1 ;
402  for (int i=0; i<nr-nl-1; i++) {
403  der_c += pari*(2*i+1)*coef(0, k, j, i) ;
404  pari = -pari ;
405  }
406  rhs.set(index) = der_c / alpha ;
407  }
408  index++ ;
409  }
410  if (base_r == R_CHEBP) {
411  for (int i=0; i<nr-nl; i++) {
412  val_b += coef(0, k, j, i) ;
413  der_b += 4*i*i*coef(0, k, j, i) ;
414  }
415  }
416  else {
417  assert(base_r == R_CHEBI) ;
418  for (int i=0; i<nr-nl-1; i++) {
419  val_b += coef(0, k, j, i) ;
420  der_b += (2*i+1)*(2*i+1)*coef(0, k, j, i) ;
421  }
422  }
423  if (domain_bc==0) {
424  rhs.set(index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha ;
425  index++ ;
426  }
427  else {
428  rhs.set(index) = -val_b ;
429  if (derivative)
430  rhs.set(index+1) = -der_b/alpha ;
431 
432  // Now for l=1
433  alpha = mp_aff->get_alpha()[1] ;
434  if ((1 == domain_bc)&&(bc_ced))
435  alpha = -0.25/alpha ;
436  nr = mgrid.get_nr(1) ;
437  int nr_lim = nr - n_conditions ;
438  val_b = 0 ; pari = 1 ;
439  for (int i=0; i<nr_lim; i++) {
440  val_b += pari*coef(1, k, j, i) ;
441  pari = -pari ;
442  }
443  rhs.set(index) += val_b ;
444  index++ ;
445  if (derivative) {
446  der_b = 0 ; pari = -1 ;
447  for (int i=0; i<nr_lim; i++) {
448  der_b += pari*i*i*coef(1, k, j, i) ;
449  pari = -pari ;
450  }
451  rhs.set(index) += der_b/alpha ;
452  index++ ;
453  }
454  val_b = 0 ;
455  for (int i=0; i<nr_lim; i++)
456  val_b += coef(1, k, j, i) ;
457  der_b = 0 ;
458  for (int i=0; i<nr_lim; i++) {
459  der_b += i*i*coef(1, k, j, i) ;
460  }
461  if (domain_bc==1) {
462  rhs.set(index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha ;
463  index++ ;
464  }
465  else {
466  rhs.set(index) = -val_b ;
467  if (derivative) rhs.set(index+1) = -der_b / alpha ;
468  }
469  }
470  }
471  else{
472  alpha = mp_aff->get_alpha()[1] ;
473  if ((1 == domain_bc)&&(bc_ced))
474  alpha = -0.25/alpha ;
475  nr = mgrid.get_nr(1) ;
476  int nr_lim = nr - 1 ;
477  val_b = 0 ;
478  for (int i=0; i<nr_lim; i++)
479  val_b += coef(1, k, j, i) ;
480  der_b = 0 ;
481  for (int i=0; i<nr_lim; i++) {
482  der_b += i*i*coef(1, k, j, i) ;
483  }
484  if (domain_bc==1) {
485  rhs.set(index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha ;
486  index++ ;
487  }
488  else {
489  rhs.set(index) = -val_b ;
490  if (derivative) rhs.set(index+1) = -der_b / alpha ;
491  }
492  }
493 
494 
495 
496  for (int lz=2; lz<=domain_bc; lz++) {
497  alpha = mp_aff->get_alpha()[lz] ;
498  if ((lz == domain_bc)&&(bc_ced))
499  alpha = -0.25/alpha ;
500  nr = mgrid.get_nr(lz) ;
501  int nr_lim = nr - n_conditions ;
502  val_b = 0 ; pari = 1 ;
503  for (int i=0; i<nr_lim; i++) {
504  val_b += pari*coef(lz, k, j, i) ;
505  pari = -pari ;
506  }
507  rhs.set(index) += val_b ;
508  index++ ;
509  if (derivative) {
510  der_b = 0 ; pari = -1 ;
511  for (int i=0; i<nr_lim; i++) {
512  der_b += pari*i*i*coef(lz, k, j, i) ;
513  pari = -pari ;
514  }
515  rhs.set(index) += der_b/alpha ;
516  index++ ;
517  }
518  val_b = 0 ;
519  for (int i=0; i<nr_lim; i++)
520  val_b += coef(lz, k, j, i) ;
521  der_b = 0 ;
522  for (int i=0; i<nr_lim; i++) {
523  der_b += i*i*coef(lz, k, j, i) ;
524  }
525  if (domain_bc==lz) {
526  rhs.set(index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha ;
527  index++ ;
528  }
529  else {
530  rhs.set(index) = -val_b ;
531  if (derivative) rhs.set(index+1) = -der_b / alpha ;
532  }
533  }
534  Tbl solut(sys_size);
535  assert(l_q>=0) ;
536  switch (l_q) {
537  case 0 :
538  solut = sys_l0.inverse(rhs) ;
539  break ;
540  case 1:
541  solut = sys_l1.inverse(rhs) ;
542  break ;
543  default:
544  solut = (l_q%2 == 0 ? sys_even.inverse(rhs) :
545  sys_odd.inverse(rhs)) ;
546  }
547  index = 0 ;
548  int dec = (base_r == R_CHEBP ? 0 : 1) ;
549  nr = mgrid.get_nr(0) ;
550  if (isexcised==false){
551  if (l_q>=2) {
552  coef.set(0, k, j, nr-2-dec) = solut(index) ;
553  index++ ;
554  }
555  coef.set(0, k, j, nr-1-dec) = solut(index) ;
556  index++ ;
557  if (base_r == R_CHEBI)
558  coef.set(0, k, j, nr-1) = 0 ;
559  if (domain_bc > 0) {
560  //Pour domaine 1
561  nr = mgrid.get_nr(1) ;
562  for (nl=1; nl<=n_conditions; nl++) {
563  int ii = n_conditions - nl + 1 ;
564  coef.set(1, k, j, nr-ii) = solut(index) ;
565  index++ ;
566  }
567  }
568  }
569  else{
570  coef.set(1,k,j,nr-1)=solut(index);
571  index++;
572  }
573  for (int lz=2; lz<=domain_bc; lz++) {
574  nr = mgrid.get_nr(lz) ;
575  for (nl=1; nl<=n_conditions; nl++) {
576  int ii = n_conditions - nl + 1 ;
577  coef.set(lz, k, j, nr-ii) = solut(index) ;
578  index++ ;
579  }
580  }
581  } //End of nullite_plm
582  else {
583  for (int lz=0; lz<=domain_bc; lz++)
584  for (int i=0; i<mgrid.get_nr(lz); i++)
585  coef.set(lz, k, j, i) = 0 ;
586  }
587  } //End of loop on j
588  } //End of loop on k
589 }
590 
591 void Scalar::match_tau_star(Param& par_bc, Param* par_mat) {
592 
593  const Map_star* mp_star = dynamic_cast<const Map_star*>(mp) ;
594  assert(mp_star != 0x0) ;
595 
596  const Mg3d& mgrid = *mp_star->get_mg() ;
597  assert(mgrid.get_type_r(0) == RARE) ;
598  assert (par_bc.get_n_tbl_mod() > 0) ;
599  Tbl mu = par_bc.get_tbl_mod(0) ;
600  if (etat == ETATZERO) {
601  if (mu.get_etat() == ETATZERO)
602  return ;
603  else
604  annule_hard() ;
605  }
606 
607  int nz = mgrid.get_nzone() ;
608  int nzm1 = nz - 1 ;
609  bool ced = (mgrid.get_type_r(nzm1) == UNSURR) ;
610  assert(par_bc.get_n_int() >= 2);
611  int domain_bc = par_bc.get_int(0) ;
612  bool bc_ced = ((ced) && (domain_bc == nzm1)) ;
613 
614  int n_conditions = par_bc.get_int(1) ;
615  assert ((n_conditions==1)||(n_conditions==2)) ;
616  bool derivative = (n_conditions == 2) ;
617  int dl = 0 ; int l_min = 0 ; int excised_i =0;
618  if (par_bc.get_n_int() > 2) {
619  dl = par_bc.get_int(2) ;
620  l_min = par_bc.get_int(3) ;
621  if(par_bc.get_n_int() >4){
622  excised_i = par_bc.get_int(4);
623  }
624  }
625  bool isexcised = (excised_i==1);
626 
627  int nt = mgrid.get_nt(0) ;
628  int np = mgrid.get_np(0) ;
629  assert (par_bc.get_n_double_mod() == 2) ;
630  double c_val = par_bc.get_double_mod(0) ;
631  double c_der = par_bc.get_double_mod(1) ;
632  if (bc_ced) {
633  c_val = 1 ;
634  c_der = 0 ;
635  mu.annule_hard() ;
636  }
637  if (mu.get_etat() == ETATZERO)
638  mu.annule_hard() ;
639  int system01_size =0;
640  int system_size =0;
641  if (isexcised ==false){
642  system01_size = 1 ;
643  system_size = 2 ;
644  }
645  else{
646  system01_size = -1;
647  system_size = -1;
648  }
649  for (int lz=1; lz<=domain_bc; lz++) {
650  system01_size += n_conditions ;
651  system_size += n_conditions ;
652  }
653  assert (etat != ETATNONDEF) ;
654 
655  bool need_matrices = true ;
656  if (par_mat != 0x0)
657  if (par_mat->get_n_matrice_mod() == 4)
658  need_matrices = false ;
659 
660  int npp2 = np ;
661  Mg3d mg_01(npp2, system01_size, system01_size, nt, SYM, SYM) ; Mg3d mgmg(npp2, system_size, system_size, nt, SYM, SYM) ;
662 
663  Valeur system_l0(mg_01) ;
664  Valeur system_l1(mg_01) ;
665  Valeur system_even(mgmg) ;
666  Valeur system_odd(mgmg) ;
667 
668 
669 
670  if (need_matrices) {
671  system_l0.annule_hard() ;
672  system_l1.annule_hard() ;
673  system_even.annule_hard() ;
674  system_odd.annule_hard() ;
675 
676  int nr = mgrid.get_nr(0);
677  Valeur alpha = mp_star->get_alpha() ;
678  for (int k=0; k<np; k++)
679  for (int j=0; j<nt; j++){
680  int index = 0 ; int index01 = 0 ;
681  if (isexcised == false){
682  system_even.set(k,j,index, index) = 1. ;
683  system_even.set(k,j,index, index + 1) = -1. ; //C_{N-1} - C_N = \sum (-1)^i C_i
684  system_odd.set(k,j,index, index) = -(2.*nr - 5.)/alpha(0,k,j,0) ;
685  system_odd.set(k,j,index, index+1) = (2.*nr - 3.)/alpha(0,k,j,0) ;
686  index++ ;
687  if (domain_bc == 0) {
688  system_l0.set(k,j,index01, index01) = c_val + c_der*4.*(nr-1)*(nr-1)/alpha(0,k,j,0) ;
689  system_l1.set(k,j,index01, index01) = c_val + c_der*(2*nr-3)*(2*nr-3)/alpha(0,k,j,0) ;
690  index01++ ;
691  system_even.set(k,j,index, index-1) = c_val + c_der*4.*(nr-2)*(nr-2)/alpha(0,k,j,0) ;
692  system_even.set(k,j,index, index) = c_val + c_der*4.*(nr-1)*(nr-1)/alpha(0,k,j,0) ;
693  system_odd.set(k,j,index, index-1) = c_val + c_der*(2*nr-5)*(2*nr-5)/alpha(0,k,j,0) ;
694  system_odd.set(k,j,index, index) = c_val + c_der*(2*nr-3)*(2*nr-3)/alpha(0,k,j,0) ;
695  index++ ;
696  }
697  else {
698  system_l0.set(k,j,index01, index01) = 1. ;
699  system_l1.set(k,j,index01, index01) = 1. ;
700  system_even.set(k,j,index, index-1) = 1. ;
701  system_even.set(k,j,index, index) = 1. ;
702  system_odd.set(k,j,index, index-1) = 1. ;
703  system_odd.set(k,j,index, index) = 1. ;
704  if (derivative) {
705  system_l0.set(k,j,index01+1, index01) = 4*(nr-1)*(nr-1)/alpha(0,k,j,0) ;
706  system_l1.set(k,j,index01+1, index01) = (2*nr-3)*(2*nr-3)/alpha(0,k,j,0) ;
707  system_even.set(k,j,index+1, index-1) = 4*(nr-2)*(nr-2)/alpha(0,k,j,0) ;
708  system_even.set(k,j,index+1, index) = 4*(nr-1)*(nr-1)/alpha(0,k,j,0) ;
709  system_odd.set(k,j,index+1, index-1) = (2*nr-5)*(2*nr-5)/alpha(0,k,j,0) ;
710  system_odd.set(k,j,index+1, index) = (2*nr-3)*(2*nr-3)/alpha(0,k,j,0) ;
711  }
712  }
713  if (domain_bc > 0) {
714  // Do things for lz=1;
715  nr = mgrid.get_nr(1) ;
716  // alpha = mp_star->get_alpha() ;
717  // if ((1 == domain_bc)&&(bc_ced))
718  // alpha = -0.25/alpha ; // No ced in Map_star
719  system_l0.set(k,j,index01, index01+1) = 1. ;
720  system_l0.set(k,j,index01, index01+2) = -1. ;
721  system_l1.set(k,j,index01, index01+1) = 1. ;
722  system_l1.set(k,j,index01, index01+2) = -1. ;
723  index01++ ;
724  system_even.set(k,j,index, index+1) = 1. ;
725  system_even.set(k,j,index, index+2) = -1. ;
726  system_odd.set(k,j,index, index+1) = 1. ;
727  system_odd.set(k,j,index, index+2) = -1. ;
728  index++ ;
729  if (derivative) {
730  system_l0.set(k,j,index01, index01) = -(nr-2)*(nr-2)/alpha(1,k,j,0) ;
731  system_l0.set(k,j,index01, index01+1) = (nr-1)*(nr-1)/alpha(1,k,j,0) ;
732  system_l1.set(k,j,index01, index01) = -(nr-2)*(nr-2)/alpha(1,k,j,0) ;
733  system_l1.set(k,j,index01, index01+1) = (nr-1)*(nr-1)/alpha(1,k,j,0) ;
734  index01++ ;
735  system_even.set(k,j,index, index) = -(nr-2)*(nr-2)/alpha(1,k,j,0) ;
736  system_even.set(k,j,index, index+1) = (nr-1)*(nr-1)/alpha(1,k,j,0) ;
737  system_odd.set(k,j,index, index) = -(nr-2)*(nr-2)/alpha(1,k,j,0) ;
738  system_odd.set(k,j,index, index+1) = (nr-1)*(nr-1)/alpha(1,k,j,0) ;
739  index++ ;
740  }
741 
742  if (1 == domain_bc) {
743  system_l0.set(k,j,index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha(1,k,j,0) ;
744  system_l0.set(k,j,index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha(1,k,j,0) ;
745  system_l1.set(k,j,index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha(1,k,j,0) ;
746  system_l1.set(k,j,index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha(1,k,j,0) ;
747  index01++ ;
748  system_even.set(k,j,index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha(1,k,j,0) ;
749  system_even.set(k,j,index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha(1,k,j,0) ;
750  system_odd.set(k,j,index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha(1,k,j,0) ;
751  system_odd.set(k,j,index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha(1,k,j,0) ;
752  index++ ;
753  }
754  else {
755  system_l0.set(k,j,index01, index01-1) = 1. ;
756  system_l0.set(k,j,index01, index01) = 1. ;
757  system_l1.set(k,j,index01, index01-1) = 1. ;
758  system_l1.set(k,j,index01, index01) = 1. ;
759  system_even.set(k,j,index, index-1) = 1. ;
760  system_even.set(k,j,index, index) = 1. ;
761  system_odd.set(k,j,index, index-1) = 1. ;
762  system_odd.set(k,j,index, index) = 1. ;
763  if (derivative) {
764  system_l0.set(k,j,index01+1, index01-1) = (nr-2)*(nr-2)/alpha(1,k,j,0) ;
765  system_l0.set(k,j,index01+1, index01) = (nr-1)*(nr-1)/alpha(1,k,j,0) ;
766  system_l1.set(k,j,index01+1, index01-1) = (nr-2)*(nr-2)/alpha(1,k,j,0) ;
767  system_l1.set(k,j,index01+1, index01) = (nr-1)*(nr-1)/alpha(1,k,j,0) ;
768  system_even.set(k,j,index+1, index-1) = (nr-2)*(nr-2)/alpha(1,k,j,0) ;
769  system_even.set(k,j,index+1, index) = (nr-1)*(nr-1)/alpha(1,k,j,0) ;
770  system_odd.set(k,j,index+1, index-1) = (nr-2)*(nr-2)/alpha(1,k,j,0) ;
771  system_odd.set(k,j,index+1, index) = (nr-1)*(nr-1)/alpha(1,k,j,0) ;
772  }
773  }
774  }
775  }
776  else {
777  nr = mgrid.get_nr(1) ;
778  alpha = mp_star->get_alpha() ;
779  // if ((1 == domain_bc)&&(bc_ced))
780  // alpha = -0.25/alpha ; // No ced in Map_star
781  if (1 == domain_bc) {
782  system_l0.set(k,j,index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha(1,k,j,0) ;
783  system_l1.set(k,j,index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha(1,k,j,0) ;
784  index01++ ;
785  system_even.set(k,j,index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha(1,k,j,0) ;
786  system_odd.set(k,j,index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha(1,k,j,0) ;
787  index++ ;
788  }
789  else {
790  system_l0.set(k,j,index01, index01) = 1. ;
791  system_l1.set(k,j,index01, index01) = 1. ;
792  system_even.set(k,j,index, index) = 1. ;
793  system_odd.set(k,j,index, index) = 1. ;
794  if (derivative) {
795  system_l0.set(k,j,index01+1, index01) = (nr-1)*(nr-1)/alpha(1,k,j,0) ;
796  system_l1.set(k,j,index01+1, index01) = (nr-1)*(nr-1)/alpha(1,k,j,0) ;
797  system_even.set(k,j,index+1, index) = (nr-1)*(nr-1)/alpha(1,k,j,0) ;
798  system_odd.set(k,j,index+1, index) = (nr-1)*(nr-1)/alpha(1,k,j,0) ;
799  }
800 
801  }
802  }
803  for (int lz=2; lz<=domain_bc; lz++) {
804  nr = mgrid.get_nr(lz) ;
805  alpha = mp_star->get_alpha() ;
806  // if ((lz == domain_bc)&&(bc_ced))
807  // alpha = -0.25/alpha ; // No ced in Map_star
808  system_l0.set(k,j,index01, index01+1) = 1. ;
809  system_l0.set(k,j,index01, index01+2) = -1. ;
810  system_l1.set(k,j,index01, index01+1) = 1. ;
811  system_l1.set(k,j,index01, index01+2) = -1. ;
812  index01++ ;
813  system_even.set(k,j,index, index+1) = 1. ;
814  system_even.set(k,j,index, index+2) = -1. ;
815  system_odd.set(k,j,index, index+1) = 1. ;
816  system_odd.set(k,j,index, index+2) = -1. ;
817  index++ ;
818  if (derivative) {
819  system_l0.set(k,j,index01, index01) = -(nr-2)*(nr-2)/alpha(lz,k,j,0) ;
820  system_l0.set(k,j,index01, index01+1) = (nr-1)*(nr-1)/alpha(lz,k,j,0) ;
821  system_l1.set(k,j,index01, index01) = -(nr-2)*(nr-2)/alpha(lz,k,j,0) ;
822  system_l1.set(k,j,index01, index01+1) = (nr-1)*(nr-1)/alpha(lz,k,j,0) ;
823  index01++ ;
824  system_even.set(k,j,index, index) = -(nr-2)*(nr-2)/alpha(lz,k,j,0) ;
825  system_even.set(k,j,index, index+1) = (nr-1)*(nr-1)/alpha(lz,k,j,0) ;
826  system_odd.set(k,j,index, index) = -(nr-2)*(nr-2)/alpha(lz,k,j,0) ;
827  system_odd.set(k,j,index, index+1) = (nr-1)*(nr-1)/alpha(lz,k,j,0) ;
828  index++ ;
829  }
830 
831  if (lz == domain_bc) {
832  system_l0.set(k,j,index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha(lz,k,j,0) ;
833  system_l0.set(k,j,index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha(lz,k,j,0) ;
834  system_l1.set(k,j,index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha(lz,k,j,0) ;
835  system_l1.set(k,j,index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha(lz,k,j,0) ;
836  index01++ ;
837  system_even.set(k,j,index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha(lz,k,j,0) ;
838  system_even.set(k,j,index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha(lz,k,j,0) ;
839  system_odd.set(k,j,index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha(lz,k,j,0) ;
840  system_odd.set(k,j,index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha(lz,k,j,0) ;
841  index++ ;
842  }
843  else {
844  system_l0.set(k,j,index01, index01-1) = 1. ;
845  system_l0.set(k,j,index01, index01) = 1. ;
846  system_l1.set(k,j,index01, index01-1) = 1. ;
847  system_l1.set(k,j,index01, index01) = 1. ;
848  system_even.set(k,j,index, index-1) = 1. ;
849  system_even.set(k,j,index, index) = 1. ;
850  system_odd.set(k,j,index, index-1) = 1. ;
851  system_odd.set(k,j,index, index) = 1. ;
852  if (derivative) {
853  system_l0.set(k,j,index01+1, index01-1) = (nr-2)*(nr-2)/alpha(lz,k,j,0) ;
854  system_l0.set(k,j,index01+1, index01) = (nr-1)*(nr-1)/alpha(lz,k,j,0) ;
855  system_l1.set(k,j,index01+1, index01-1) = (nr-2)*(nr-2)/alpha(lz,k,j,0) ;
856  system_l1.set(k,j,index01+1, index01) = (nr-1)*(nr-1)/alpha(lz,k,j,0) ;
857  system_even.set(k,j,index+1, index-1) = (nr-2)*(nr-2)/alpha(lz,k,j,0) ;
858  system_even.set(k,j,index+1, index) = (nr-1)*(nr-1)/alpha(lz,k,j,0) ;
859  system_odd.set(k,j,index+1, index-1) = (nr-2)*(nr-2)/alpha(lz,k,j,0) ;
860  system_odd.set(k,j,index+1, index) = (nr-1)*(nr-1)/alpha(lz,k,j,0) ;
861  }
862  }
863  } // End of loop on lz
864 
865  assert(index01 == system01_size) ;
866  assert(index == system_size) ;
867  // system_l0.set_lu() ;
868  // system_l1.set_lu() ;
869  // system_even.set_lu() ;
870  // system_odd.set_lu() ;
871  // if (par_mat != 0x0) {
872  // Matrice* pmat0 = new Matrice(system_l0) ;
873  // Matrice* pmat1 = new Matrice(system_l1) ;
874  // Matrice* pmat2 = new Matrice(system_even) ;
875  // Matrice* pmat3 = new Matrice(system_odd) ;
876  // par_mat->add_matrice_mod(*pmat0, 0) ;
877  // par_mat->add_matrice_mod(*pmat1, 1) ;
878  // par_mat->add_matrice_mod(*pmat2, 2) ;
879  // par_mat->add_matrice_mod(*pmat3, 3) ;
880  // }
881  }
882  }// End of if (need_matrices) ...
883 
884  const Valeur& sys_l0 = system_l0 ; //(need_matrices ? system_l0 : par_mat->get_matrice_mod(0) ) ;
885  const Valeur& sys_l1 = system_l1 ; //(need_matrices ? system_l1 : par_mat->get_matrice_mod(1) ) ;
886  const Valeur& sys_even = system_even ; //(need_matrices ? system_even :
887  //par_mat->get_matrice_mod(2) ) ;
888  const Valeur& sys_odd = system_odd ; //(need_matrices ? system_odd :
889  //par_mat->get_matrice_mod(3) ) ;
890 
891  va.coef() ;
892  va.ylm() ;
893  const Base_val& base = get_spectral_base() ;
894  Mtbl_cf& coef = *va.c_cf ; cout << "COEF : " << coef << endl;
895  if (va.c != 0x0) {
896  delete va.c ;
897  va.c = 0x0 ;
898  }
899  int m_q, l_q, base_r ;
900  for (int k=0; k<np+2; k++) {
901  for (int j=0; j<nt; j++) {
902  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;//#0 here as domain index
903  if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q >= l_min)) {
904  l_q += dl ;
905  int sys_size = (l_q < 2 ? system01_size : system_size) ;
906  int nl = (l_q < 2 ? 1 : 2) ;
907  Tbl rhs(np+2,nt,sys_size) ; rhs.annule_hard() ;
908  int index = 0 ;
909  int pari = 1 ;
910  Valeur alpha = mp_star->get_alpha() ;
911  int nr = mgrid.get_nr(0) ;
912  double val_b = 0. ;
913  double der_b =0. ;
914  if (isexcised==false){
915  if (l_q>=2) {
916  if (base_r == R_CHEBP) {
917  double val_c = 0.; pari = 1 ;
918  for (int i=0; i<nr-nl; i++) {
919  val_c += pari*coef(0, k, j, i) ;
920  pari = -pari ;
921  }
922  rhs.set(k,j,index) = val_c ;
923  }
924  else {
925  assert( base_r == R_CHEBI) ;
926  double der_c = 0.; pari = 1 ;
927  for (int i=0; i<nr-nl-1; i++) {
928  der_c += pari*(2*i+1)*coef(0, k, j, i) ;
929  pari = -pari ;
930  }
931  rhs.set(k,j,index) = der_c / alpha(0, k, j, 0) ;
932  }
933  index++ ;
934  }
935  if (base_r == R_CHEBP) {
936  for (int i=0; i<nr-nl; i++) {
937  val_b += coef(0, k, j, i) ;
938  der_b += 4*i*i*coef(0, k, j, i) ;
939  }
940  }
941  else {
942  assert(base_r == R_CHEBI) ;
943  for (int i=0; i<nr-nl-1; i++) {
944  val_b += coef(0, k, j, i) ;
945  der_b += (2*i+1)*(2*i+1)*coef(0, k, j, i) ;
946  }
947  }
948  if (domain_bc==0) {
949  rhs.set(k,j,index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha(0, k, j, 0) ;
950  index++ ;
951  }
952  else {
953  rhs.set(k,j,index) = -val_b ;
954  if (derivative)
955  rhs.set(k,j,index+1) = -der_b/alpha(0, k, j, 0) ;
956 
957  // Now for l=1
958  alpha = mp_star->get_alpha() ;
959  // if ((1 == domain_bc)&&(bc_ced))
960  // alpha = -0.25/alpha ; // No ced in Map_star
961  nr = mgrid.get_nr(1) ;
962  int nr_lim = nr - n_conditions ;
963  val_b = 0 ; pari = 1 ;
964  for (int i=0; i<nr_lim; i++) {
965  val_b += pari*coef(1, k, j, i) ;
966  pari = -pari ;
967  }
968  rhs.set(k,j,index) += val_b ;
969  index++ ;
970  if (derivative) {
971  der_b = 0 ; pari = -1 ;
972  for (int i=0; i<nr_lim; i++) {
973  der_b += pari*i*i*coef(1, k, j, i) ;
974  pari = -pari ;
975  }
976  rhs.set(k,j,index) += der_b/alpha(1,k,j,0) ;
977  index++ ;
978  }
979  val_b = 0 ;
980  for (int i=0; i<nr_lim; i++)
981  val_b += coef(1, k, j, i) ;
982  der_b = 0 ;
983  for (int i=0; i<nr_lim; i++) {
984  der_b += i*i*coef(1, k, j, i) ;
985  }
986  if (domain_bc==1) {
987  rhs.set(k,j,index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha(1,k,j,0) ;
988  index++ ;
989  }
990  else {
991  rhs.set(k,j,index) = -val_b ;
992  if (derivative) rhs.set(k,j,index+1) = -der_b / alpha(1,k,j,0) ;
993  }
994  }
995  }
996  else{
997  alpha = mp_star->get_alpha() ;
998  // if ((1 == domain_bc)&&(bc_ced))
999  // alpha = -0.25/alpha ; // No ced in Map_star
1000  nr = mgrid.get_nr(1) ;
1001  int nr_lim = nr - 1 ;
1002  val_b = 0 ;
1003  for (int i=0; i<nr_lim; i++)
1004  val_b += coef(1, k, j, i) ;
1005  der_b = 0 ;
1006  for (int i=0; i<nr_lim; i++) {
1007  der_b += i*i*coef(1, k, j, i) ;
1008  }
1009  if (domain_bc==1) {
1010  rhs.set(k,j,index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha(1,k,j,0) ;
1011  index++ ;
1012  }
1013  else {
1014  rhs.set(k,j,index) = -val_b ;
1015  if (derivative) rhs.set(k,j,index+1) = -der_b / alpha(1,k,j,0) ;
1016  }
1017  }
1018 
1019 
1020 
1021  for (int lz=2; lz<=domain_bc; lz++) {
1022  alpha = mp_star->get_alpha() ;
1023  // if ((lz == domain_bc)&&(bc_ced))
1024  // alpha = -0.25/alpha ; No ced in Map_star
1025  nr = mgrid.get_nr(lz) ;
1026  int nr_lim = nr - n_conditions ;
1027  val_b = 0 ; pari = 1 ;
1028  for (int i=0; i<nr_lim; i++) {
1029  val_b += pari*coef(lz, k, j, i) ;
1030  pari = -pari ;
1031  }
1032  rhs.set(k,j,index) += val_b ;
1033  index++ ;
1034  if (derivative) {
1035  der_b = 0 ; pari = -1 ;
1036  for (int i=0; i<nr_lim; i++) {
1037  der_b += pari*i*i*coef(lz, k, j, i) ;
1038  pari = -pari ;
1039  }
1040  rhs.set(k,j,index) += der_b/alpha(lz, k, j, 0) ;
1041  index++ ;
1042  }
1043  val_b = 0 ;
1044  for (int i=0; i<nr_lim; i++)
1045  val_b += coef(lz, k, j, i) ;
1046  der_b = 0 ;
1047  for (int i=0; i<nr_lim; i++) {
1048  der_b += i*i*coef(lz, k, j, i) ;
1049  }
1050  if (domain_bc==lz) {
1051  rhs.set(k,j,index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha(lz, k, j, 0) ;
1052  index++ ;
1053  }
1054  else {
1055  rhs.set(k,j,index) = -val_b ;
1056  if (derivative) rhs.set(k,j,index+1) = -der_b / alpha(lz, k, j, 0) ;
1057  }
1058  }
1059  Tbl solut(sys_size);
1060  assert(l_q>=0) ;
1061  Tbl rhs_tmp(sys_size) ; rhs_tmp.set_etat_qcq() ;
1062  switch (l_q) {
1063  case 0 :{
1064  Matrice sys_l0_tmp(sys_size, sys_size) ; sys_l0_tmp.set_etat_qcq() ;
1065  for (int ind1=0; ind1<sys_size; ind1++){
1066  rhs_tmp.set(ind1) = rhs(k,j,ind1) ;
1067  for (int ind2=0; ind2<sys_size; ind2++)
1068  sys_l0_tmp.set(ind1,ind2) = sys_l0(k,j,ind1,ind2) ;
1069  }
1070  sys_l0_tmp.set_lu() ;
1071  solut = sys_l0_tmp.inverse(rhs_tmp) ;
1072  }
1073  break ;
1074  case 1:{
1075  Matrice sys_l1_tmp(sys_size, sys_size) ; sys_l1_tmp.set_etat_qcq() ;
1076  for (int ind1=0; ind1<sys_size; ind1++){
1077  rhs_tmp.set(ind1) = rhs(k,j,ind1) ;
1078  for (int ind2=0; ind2<sys_size; ind2++)
1079  sys_l1_tmp.set(ind1,ind2) = sys_l1(k,j,ind1,ind2) ;
1080  }
1081  sys_l1_tmp.set_lu() ;
1082  solut = sys_l1_tmp.inverse(rhs_tmp) ;
1083  }
1084  break ;
1085  default:
1086  {
1087  if (l_q%2 == 0){
1088  Matrice sys_even_tmp(sys_size, sys_size) ; sys_even_tmp.set_etat_qcq() ;
1089  for (int ind1=0; ind1<sys_size; ind1++){
1090  rhs_tmp.set(ind1) = rhs(k,j,ind1) ;
1091  for (int ind2=0; ind2<sys_size; ind2++)
1092  sys_even_tmp.set(ind1,ind2) = sys_even(k,j,ind1,ind2) ;
1093  }
1094  sys_even_tmp.set_lu() ;
1095  solut = sys_even_tmp.inverse(rhs_tmp) ;
1096  }
1097  else{
1098  Matrice sys_odd_tmp(sys_size, sys_size) ; sys_odd_tmp.set_etat_qcq() ;
1099  for (int ind1=0; ind1<sys_size; ind1++){
1100  rhs_tmp.set(ind1) = rhs(k,j,ind1) ;
1101  for (int ind2=0; ind2<sys_size; ind2++)
1102  sys_odd_tmp.set(ind1,ind2) = sys_odd(k,j,ind1,ind2) ;
1103  }
1104  sys_odd_tmp.set_lu() ;
1105  solut = sys_odd_tmp.inverse(rhs_tmp) ;
1106  }
1107  }
1108  }
1109  index = 0 ;
1110  int dec = (base_r == R_CHEBP ? 0 : 1) ;
1111  nr = mgrid.get_nr(0) ;
1112  if (isexcised==false){
1113  if (l_q>=2) {
1114  coef.set(0, k, j, nr-2-dec) = solut(index) ;
1115  index++ ;
1116  }
1117  coef.set(0, k, j, nr-1-dec) = solut(index) ;
1118  index++ ;
1119  if (base_r == R_CHEBI)
1120  coef.set(0, k, j, nr-1) = 0 ;
1121  if (domain_bc > 0) {
1122  //Pour domaine 1
1123  nr = mgrid.get_nr(1) ;
1124  for (nl=1; nl<=n_conditions; nl++) {
1125  int ii = n_conditions - nl + 1 ;
1126  coef.set(1, k, j, nr-ii) = solut(index) ;
1127  index++ ;
1128  }
1129  }
1130  }
1131  else{
1132  coef.set(1,k,j,nr-1)=solut(index);
1133  index++;
1134  }
1135  for (int lz=2; lz<=domain_bc; lz++) {
1136  nr = mgrid.get_nr(lz) ;
1137  for (nl=1; nl<=n_conditions; nl++) {
1138  int ii = n_conditions - nl + 1 ;
1139  coef.set(lz, k, j, nr-ii) = solut(index) ;
1140  index++ ;
1141  }
1142  }
1143  } //End of nullite_plm
1144  else {
1145  for (int lz=0; lz<=domain_bc; lz++)
1146  for (int i=0; i<mgrid.get_nr(lz); i++)
1147  coef.set(lz, k, j, i) = 0 ;
1148  }
1149  } //End of loop on j
1150  } //End of loop on k
1151 }
1152 
1153 }
Lorene::Scalar::get_spectral_base
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Definition: scalar.h:1328
Lorene::Map_star
Definition: map_star.h:23
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
Lorene::Param::get_double_mod
double & get_double_mod(int position=0) const
Returns the reference of a stored modifiable double .
Definition: param.C:501
Lorene::Map_af::get_alpha
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:607
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:151
Lorene::Param::get_n_matrice_mod
int get_n_matrice_mod() const
Returns the number of modifiable Matrice &#39;s addresses in the list.
Definition: param.C:862
Lorene::Base_val::give_quant_numbers
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Definition: base_val_quantum.C:84
Lorene::Matrice::set_lu
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition: matrice.C:395
Lorene
Lorene prototypes.
Definition: app_hor.h:67
Lorene::Mtbl_cf::set
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:304
Lorene::Matrice::inverse
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition: matrice.C:427
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:141
Lorene::Valeur::annule_hard
void annule_hard()
Sets the Valeur to zero in a hard way.
Definition: valeur.C:726
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Lorene::Param::add_matrice_mod
void add_matrice_mod(Matrice &ti, int position=0)
Adds the address of a new modifiable Matrice to the list.
Definition: param.C:869
Lorene::Scalar::match_tau
void match_tau(Param &par_bc, Param *par_mat=0x0)
Method for matching accross domains and imposing boundary condition.
Definition: scalar_match_tau.C:75
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
Lorene::Scalar::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:386
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
Lorene::Param::get_n_int
int get_n_int() const
Returns the number of stored int &#39;s addresses.
Definition: param.C:242
Lorene::Param::get_int
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:295
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Lorene::Matrice
Matrix handling.
Definition: matrice.h:152
Lorene::Valeur::c
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:309
Lorene::Param
Parameter storage.
Definition: param.h:125
Lorene::Param::get_tbl_mod
Tbl & get_tbl_mod(int position=0) const
Returns the reference of a modifiable Tbl stored in the list.
Definition: param.C:639
Lorene::Param::get_n_double_mod
int get_n_double_mod() const
Returns the number of stored modifiable double &#39;s addresses.
Definition: param.C:449
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Lorene::Scalar::va
Valeur va
The numerical value of the Scalar.
Definition: scalar.h:411
Lorene::Matrice::annule_hard
void annule_hard()
Sets the logical state to ETATQCQ (undefined state).
Definition: matrice.C:196
Lorene::Map_star::get_alpha
const Valeur & get_alpha() const
Returns the reference on the Tbl alpha.
Definition: map_star.C:237
Lorene::Matrice::set
double & set(int j, int i)
Read/write of a particuliar element.
Definition: matrice.h:277
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Lorene::Param::get_n_tbl_mod
int get_n_tbl_mod() const
Returns the number of modifiable Tbl &#39;s addresses in the list.
Definition: param.C:587
Lorene::Mg3d
Multi-domain grid.
Definition: grilles.h:279
Lorene::Base_val
Bases of the spectral expansions.
Definition: base_val.h:325
Lorene::Map_af
Affine radial mapping.
Definition: map.h:2048
Lorene::Scalar::etat
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary)...
Definition: scalar.h:402
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:196
Lorene::Matrice::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:178
Lorene::Param::get_matrice_mod
Matrice & get_matrice_mod(int position=0) const
Returns the reference of a modifiable Matrice stored in the list.
Definition: param.C:914
Lorene::Tbl
Basic array class.
Definition: tbl.h:164
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Lorene::Mg3d::get_type_r
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:491
Lorene::Scalar::match_tau_star
void match_tau_star(Param &par_bc, Param *par_mat=0x0)
Method for matching accross domains and imposing boundary condition.
Definition: scalar_match_tau.C:591
Lorene::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:301
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:375
Lorene::Valeur::set
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition: valeur.h:373