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