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.7 2016/12/05 16:18:18 j_novak Exp $
31  * $Log: scalar_match_tau.C,v $
32  * Revision 1.7 2016/12/05 16:18:18 j_novak
33  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
34  *
35  * Revision 1.6 2014/10/13 08:53:46 j_novak
36  * Lorene classes and functions now belong to the namespace Lorene.
37  *
38  * Revision 1.5 2014/10/06 15:16:16 j_novak
39  * Modified #include directives to use c++ syntax.
40  *
41  * Revision 1.4 2012/05/11 14:11:24 j_novak
42  * Modifications to deal with the accretion of a scalar field
43  *
44  * Revision 1.3 2012/02/06 12:59:07 j_novak
45  * Correction of some errors.
46  *
47  * Revision 1.2 2008/08/20 13:23:43 j_novak
48  * The shift in quantum number l (e.g. for \tilde{B}) is now taken into account.
49  *
50  * Revision 1.1 2008/05/24 15:05:22 j_novak
51  * New method Scalar::match_tau to match the output of an explicit time-marching scheme with the tau method.
52  *
53  *
54  * $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_match_tau.C,v 1.7 2016/12/05 16:18:18 j_novak Exp $
55  *
56  */
57 
58 // C headers
59 #include <cassert>
60 #include <cmath>
61 
62 // Lorene headers
63 #include "tensor.h"
64 #include "matrice.h"
65 #include "proto.h"
66 #include "param.h"
67 
68 namespace Lorene {
69 void Scalar::match_tau(Param& par_bc, Param* par_mat) {
70 
71  const Map_af* mp_aff = dynamic_cast<const Map_af*>(mp) ;
72  assert(mp_aff != 0x0) ; //Only affine mapping for the moment
73 
74  const Mg3d& mgrid = *mp_aff->get_mg() ;
75  assert(mgrid.get_type_r(0) == RARE) ;
76  assert (par_bc.get_n_tbl_mod() > 0) ;
77  Tbl mu = par_bc.get_tbl_mod(0) ;
78  if (etat == ETATZERO) {
79  if (mu.get_etat() == ETATZERO)
80  return ;
81  else
82  annule_hard() ;
83  }
84 
85  int nz = mgrid.get_nzone() ;
86  int nzm1 = nz - 1 ;
87  bool ced = (mgrid.get_type_r(nzm1) == UNSURR) ;
88  assert(par_bc.get_n_int() >= 2);
89  int domain_bc = par_bc.get_int(0) ;
90  bool bc_ced = ((ced) && (domain_bc == nzm1)) ;
91 
92  int n_conditions = par_bc.get_int(1) ;
93  assert ((n_conditions==1)||(n_conditions==2)) ;
94  bool derivative = (n_conditions == 2) ;
95  int dl = 0 ; int l_min = 0 ; int excised_i =0;
96  if (par_bc.get_n_int() > 2) {
97  dl = par_bc.get_int(2) ;
98  l_min = par_bc.get_int(3) ;
99  if(par_bc.get_n_int() >4){
100  excised_i = par_bc.get_int(4);
101  }
102  }
103  bool isexcised = (excised_i==1);
104 
105  int nt = mgrid.get_nt(0) ;
106  int np = mgrid.get_np(0) ;
107  assert (par_bc.get_n_double_mod() == 2) ;
108  double c_val = par_bc.get_double_mod(0) ;
109  double c_der = par_bc.get_double_mod(1) ;
110  if (bc_ced) {
111  c_val = 1 ;
112  c_der = 0 ;
113  mu.annule_hard() ;
114  }
115  if (mu.get_etat() == ETATZERO)
116  mu.annule_hard() ;
117  int system01_size =0;
118  int system_size =0;
119  if (isexcised ==false){
120  system01_size = 1 ;
121  system_size = 2 ;
122  }
123  else{
124  system01_size = -1;
125  system_size = -1;
126  }
127  for (int lz=1; lz<=domain_bc; lz++) {
128  system01_size += n_conditions ;
129  system_size += n_conditions ;
130  }
131  assert (etat != ETATNONDEF) ;
132 
133  bool need_matrices = true ;
134  if (par_mat != 0x0)
135  if (par_mat->get_n_matrice_mod() == 4)
136  need_matrices = false ;
137 
138  Matrice system_l0(system01_size, system01_size) ;
139  Matrice system_l1(system01_size, system01_size) ;
140  Matrice system_even(system_size, system_size) ;
141  Matrice system_odd(system_size, system_size) ;
142 
143  if (need_matrices) {
144  system_l0.annule_hard() ;
145  system_l1.annule_hard() ;
146  system_even.annule_hard() ;
147  system_odd.annule_hard() ;
148  int index = 0 ; int index01 = 0 ;
149  int nr = mgrid.get_nr(0);
150  double alpha = mp_aff->get_alpha()[0] ;
151  if (isexcised == false){
152  system_even.set(index, index) = 1. ;
153  system_even.set(index, index + 1) = -1. ; //C_{N-1} - C_N = \sum (-1)^i C_i
154  system_odd.set(index, index) = -(2.*nr - 5.)/alpha ;
155  system_odd.set(index, index+1) = (2.*nr - 3.)/alpha ;
156  index++ ;
157  if (domain_bc == 0) {
158  system_l0.set(index01, index01) = c_val + c_der*4.*(nr-1)*(nr-1)/alpha ;
159  system_l1.set(index01, index01) = c_val + c_der*(2*nr-3)*(2*nr-3)/alpha ;
160  index01++ ;
161  system_even.set(index, index-1) = c_val + c_der*4.*(nr-2)*(nr-2)/alpha ;
162  system_even.set(index, index) = c_val + c_der*4.*(nr-1)*(nr-1)/alpha ;
163  system_odd.set(index, index-1) = c_val + c_der*(2*nr-5)*(2*nr-5)/alpha ;
164  system_odd.set(index, index) = c_val + c_der*(2*nr-3)*(2*nr-3)/alpha ;
165  index++ ;
166  }
167  else {
168  system_l0.set(index01, index01) = 1. ;
169  system_l1.set(index01, index01) = 1. ;
170  system_even.set(index, index-1) = 1. ;
171  system_even.set(index, index) = 1. ;
172  system_odd.set(index, index-1) = 1. ;
173  system_odd.set(index, index) = 1. ;
174  if (derivative) {
175  system_l0.set(index01+1, index01) = 4*(nr-1)*(nr-1)/alpha ;
176  system_l1.set(index01+1, index01) = (2*nr-3)*(2*nr-3)/alpha ;
177  system_even.set(index+1, index-1) = 4*(nr-2)*(nr-2)/alpha ;
178  system_even.set(index+1, index) = 4*(nr-1)*(nr-1)/alpha ;
179  system_odd.set(index+1, index-1) = (2*nr-5)*(2*nr-5)/alpha ;
180  system_odd.set(index+1, index) = (2*nr-3)*(2*nr-3)/alpha ;
181  }
182  }
183  if (domain_bc > 0) {
184  // Do things for lz=1;
185  nr = mgrid.get_nr(1) ;
186  alpha = mp_aff->get_alpha()[1] ;
187  if ((1 == domain_bc)&&(bc_ced))
188  alpha = -0.25/alpha ;
189  system_l0.set(index01, index01+1) = 1. ;
190  system_l0.set(index01, index01+2) = -1. ;
191  system_l1.set(index01, index01+1) = 1. ;
192  system_l1.set(index01, index01+2) = -1. ;
193  index01++ ;
194  system_even.set(index, index+1) = 1. ;
195  system_even.set(index, index+2) = -1. ;
196  system_odd.set(index, index+1) = 1. ;
197  system_odd.set(index, index+2) = -1. ;
198  index++ ;
199  if (derivative) {
200  system_l0.set(index01, index01) = -(nr-2)*(nr-2)/alpha ;
201  system_l0.set(index01, index01+1) = (nr-1)*(nr-1)/alpha ;
202  system_l1.set(index01, index01) = -(nr-2)*(nr-2)/alpha ;
203  system_l1.set(index01, index01+1) = (nr-1)*(nr-1)/alpha ;
204  index01++ ;
205  system_even.set(index, index) = -(nr-2)*(nr-2)/alpha ;
206  system_even.set(index, index+1) = (nr-1)*(nr-1)/alpha ;
207  system_odd.set(index, index) = -(nr-2)*(nr-2)/alpha ;
208  system_odd.set(index, index+1) = (nr-1)*(nr-1)/alpha ;
209  index++ ;
210  }
211 
212  if (1 == domain_bc) {
213  system_l0.set(index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
214  system_l0.set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
215  system_l1.set(index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
216  system_l1.set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
217  index01++ ;
218  system_even.set(index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
219  system_even.set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
220  system_odd.set(index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
221  system_odd.set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
222  index++ ;
223  }
224  else {
225  system_l0.set(index01, index01-1) = 1. ;
226  system_l0.set(index01, index01) = 1. ;
227  system_l1.set(index01, index01-1) = 1. ;
228  system_l1.set(index01, index01) = 1. ;
229  system_even.set(index, index-1) = 1. ;
230  system_even.set(index, index) = 1. ;
231  system_odd.set(index, index-1) = 1. ;
232  system_odd.set(index, index) = 1. ;
233  if (derivative) {
234  system_l0.set(index01+1, index01-1) = (nr-2)*(nr-2)/alpha ;
235  system_l0.set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
236  system_l1.set(index01+1, index01-1) = (nr-2)*(nr-2)/alpha ;
237  system_l1.set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
238  system_even.set(index+1, index-1) = (nr-2)*(nr-2)/alpha ;
239  system_even.set(index+1, index) = (nr-1)*(nr-1)/alpha ;
240  system_odd.set(index+1, index-1) = (nr-2)*(nr-2)/alpha ;
241  system_odd.set(index+1, index) = (nr-1)*(nr-1)/alpha ;
242  }
243  }
244  }
245  }
246  else {
247  nr = mgrid.get_nr(1) ;
248  alpha = mp_aff->get_alpha()[1] ;
249  if ((1 == domain_bc)&&(bc_ced))
250  alpha = -0.25/alpha ;
251  if (1 == domain_bc) {
252  system_l0.set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
253  system_l1.set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
254  index01++ ;
255  system_even.set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
256  system_odd.set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
257  index++ ;
258  }
259  else {
260  system_l0.set(index01, index01) = 1. ;
261  system_l1.set(index01, index01) = 1. ;
262  system_even.set(index, index) = 1. ;
263  system_odd.set(index, index) = 1. ;
264  if (derivative) {
265  system_l0.set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
266  system_l1.set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
267  system_even.set(index+1, index) = (nr-1)*(nr-1)/alpha ;
268  system_odd.set(index+1, index) = (nr-1)*(nr-1)/alpha ;
269  }
270 
271  }
272  }
273  for (int lz=2; lz<=domain_bc; lz++) {
274  nr = mgrid.get_nr(lz) ;
275  alpha = mp_aff->get_alpha()[lz] ;
276  if ((lz == domain_bc)&&(bc_ced))
277  alpha = -0.25/alpha ;
278  system_l0.set(index01, index01+1) = 1. ;
279  system_l0.set(index01, index01+2) = -1. ;
280  system_l1.set(index01, index01+1) = 1. ;
281  system_l1.set(index01, index01+2) = -1. ;
282  index01++ ;
283  system_even.set(index, index+1) = 1. ;
284  system_even.set(index, index+2) = -1. ;
285  system_odd.set(index, index+1) = 1. ;
286  system_odd.set(index, index+2) = -1. ;
287  index++ ;
288  if (derivative) {
289  system_l0.set(index01, index01) = -(nr-2)*(nr-2)/alpha ;
290  system_l0.set(index01, index01+1) = (nr-1)*(nr-1)/alpha ;
291  system_l1.set(index01, index01) = -(nr-2)*(nr-2)/alpha ;
292  system_l1.set(index01, index01+1) = (nr-1)*(nr-1)/alpha ;
293  index01++ ;
294  system_even.set(index, index) = -(nr-2)*(nr-2)/alpha ;
295  system_even.set(index, index+1) = (nr-1)*(nr-1)/alpha ;
296  system_odd.set(index, index) = -(nr-2)*(nr-2)/alpha ;
297  system_odd.set(index, index+1) = (nr-1)*(nr-1)/alpha ;
298  index++ ;
299  }
300 
301  if (lz == domain_bc) {
302  system_l0.set(index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
303  system_l0.set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
304  system_l1.set(index01, index01-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
305  system_l1.set(index01, index01) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
306  index01++ ;
307  system_even.set(index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
308  system_even.set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
309  system_odd.set(index, index-1) = c_val + c_der*(nr-2)*(nr-2)/alpha ;
310  system_odd.set(index, index) = c_val + c_der*(nr-1)*(nr-1)/alpha ;
311  index++ ;
312  }
313  else {
314  system_l0.set(index01, index01-1) = 1. ;
315  system_l0.set(index01, index01) = 1. ;
316  system_l1.set(index01, index01-1) = 1. ;
317  system_l1.set(index01, index01) = 1. ;
318  system_even.set(index, index-1) = 1. ;
319  system_even.set(index, index) = 1. ;
320  system_odd.set(index, index-1) = 1. ;
321  system_odd.set(index, index) = 1. ;
322  if (derivative) {
323  system_l0.set(index01+1, index01-1) = (nr-2)*(nr-2)/alpha ;
324  system_l0.set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
325  system_l1.set(index01+1, index01-1) = (nr-2)*(nr-2)/alpha ;
326  system_l1.set(index01+1, index01) = (nr-1)*(nr-1)/alpha ;
327  system_even.set(index+1, index-1) = (nr-2)*(nr-2)/alpha ;
328  system_even.set(index+1, index) = (nr-1)*(nr-1)/alpha ;
329  system_odd.set(index+1, index-1) = (nr-2)*(nr-2)/alpha ;
330  system_odd.set(index+1, index) = (nr-1)*(nr-1)/alpha ;
331  }
332  }
333  } // End of loop on lz
334 
335  assert(index01 == system01_size) ;
336  assert(index == system_size) ;
337  system_l0.set_lu() ;
338  system_l1.set_lu() ;
339  system_even.set_lu() ;
340  system_odd.set_lu() ;
341  if (par_mat != 0x0) {
342  Matrice* pmat0 = new Matrice(system_l0) ;
343  Matrice* pmat1 = new Matrice(system_l1) ;
344  Matrice* pmat2 = new Matrice(system_even) ;
345  Matrice* pmat3 = new Matrice(system_odd) ;
346  par_mat->add_matrice_mod(*pmat0, 0) ;
347  par_mat->add_matrice_mod(*pmat1, 1) ;
348  par_mat->add_matrice_mod(*pmat2, 2) ;
349  par_mat->add_matrice_mod(*pmat3, 3) ;
350  }
351  }// End of if (need_matrices) ...
352 
353  const Matrice& sys_l0 = (need_matrices ? system_l0 : par_mat->get_matrice_mod(0) ) ;
354  const Matrice& sys_l1 = (need_matrices ? system_l1 : par_mat->get_matrice_mod(1) ) ;
355  const Matrice& sys_even = (need_matrices ? system_even :
356  par_mat->get_matrice_mod(2) ) ;
357  const Matrice& sys_odd = (need_matrices ? system_odd :
358  par_mat->get_matrice_mod(3) ) ;
359 
360  va.coef() ;
361  va.ylm() ;
362  const Base_val& base = get_spectral_base() ;
363  Mtbl_cf& coef = *va.c_cf ;
364  if (va.c != 0x0) {
365  delete va.c ;
366  va.c = 0x0 ;
367  }
368  int m_q, l_q, base_r ;
369  for (int k=0; k<np+2; k++) {
370  for (int j=0; j<nt; j++) {
371  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;//#0 here as domain index
372  if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q >= l_min)) {
373  l_q += dl ;
374  int sys_size = (l_q < 2 ? system01_size : system_size) ;
375  int nl = (l_q < 2 ? 1 : 2) ;
376  Tbl rhs(sys_size) ; rhs.annule_hard() ;
377  int index = 0 ;
378  int pari = 1 ;
379  double alpha = mp_aff->get_alpha()[0] ;
380  int nr = mgrid.get_nr(0) ;
381  double val_b = 0. ;
382  double der_b =0. ;
383  if (isexcised==false){
384  if (l_q>=2) {
385  if (base_r == R_CHEBP) {
386  double val_c = 0.; pari = 1 ;
387  for (int i=0; i<nr-nl; i++) {
388  val_c += pari*coef(0, k, j, i) ;
389  pari = -pari ;
390  }
391  rhs.set(index) = val_c ;
392  }
393  else {
394  assert( base_r == R_CHEBI) ;
395  double der_c = 0.; pari = 1 ;
396  for (int i=0; i<nr-nl-1; i++) {
397  der_c += pari*(2*i+1)*coef(0, k, j, i) ;
398  pari = -pari ;
399  }
400  rhs.set(index) = der_c / alpha ;
401  }
402  index++ ;
403  }
404  if (base_r == R_CHEBP) {
405  for (int i=0; i<nr-nl; i++) {
406  val_b += coef(0, k, j, i) ;
407  der_b += 4*i*i*coef(0, k, j, i) ;
408  }
409  }
410  else {
411  assert(base_r == R_CHEBI) ;
412  for (int i=0; i<nr-nl-1; i++) {
413  val_b += coef(0, k, j, i) ;
414  der_b += (2*i+1)*(2*i+1)*coef(0, k, j, i) ;
415  }
416  }
417  if (domain_bc==0) {
418  rhs.set(index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha ;
419  index++ ;
420  }
421  else {
422  rhs.set(index) = -val_b ;
423  if (derivative)
424  rhs.set(index+1) = -der_b/alpha ;
425 
426  // Now for l=1
427  alpha = mp_aff->get_alpha()[1] ;
428  if ((1 == domain_bc)&&(bc_ced))
429  alpha = -0.25/alpha ;
430  nr = mgrid.get_nr(1) ;
431  int nr_lim = nr - n_conditions ;
432  val_b = 0 ; pari = 1 ;
433  for (int i=0; i<nr_lim; i++) {
434  val_b += pari*coef(1, k, j, i) ;
435  pari = -pari ;
436  }
437  rhs.set(index) += val_b ;
438  index++ ;
439  if (derivative) {
440  der_b = 0 ; pari = -1 ;
441  for (int i=0; i<nr_lim; i++) {
442  der_b += pari*i*i*coef(1, k, j, i) ;
443  pari = -pari ;
444  }
445  rhs.set(index) += der_b/alpha ;
446  index++ ;
447  }
448  val_b = 0 ;
449  for (int i=0; i<nr_lim; i++)
450  val_b += coef(1, k, j, i) ;
451  der_b = 0 ;
452  for (int i=0; i<nr_lim; i++) {
453  der_b += i*i*coef(1, k, j, i) ;
454  }
455  if (domain_bc==1) {
456  rhs.set(index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha ;
457  index++ ;
458  }
459  else {
460  rhs.set(index) = -val_b ;
461  if (derivative) rhs.set(index+1) = -der_b / alpha ;
462  }
463  }
464  }
465  else{
466  alpha = mp_aff->get_alpha()[1] ;
467  if ((1 == domain_bc)&&(bc_ced))
468  alpha = -0.25/alpha ;
469  nr = mgrid.get_nr(1) ;
470  int nr_lim = nr - 1 ;
471  val_b = 0 ;
472  for (int i=0; i<nr_lim; i++)
473  val_b += coef(1, k, j, i) ;
474  der_b = 0 ;
475  for (int i=0; i<nr_lim; i++) {
476  der_b += i*i*coef(1, k, j, i) ;
477  }
478  if (domain_bc==1) {
479  rhs.set(index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha ;
480  index++ ;
481  }
482  else {
483  rhs.set(index) = -val_b ;
484  if (derivative) rhs.set(index+1) = -der_b / alpha ;
485  }
486  }
487 
488 
489 
490  for (int lz=2; lz<=domain_bc; lz++) {
491  alpha = mp_aff->get_alpha()[lz] ;
492  if ((lz == domain_bc)&&(bc_ced))
493  alpha = -0.25/alpha ;
494  nr = mgrid.get_nr(lz) ;
495  int nr_lim = nr - n_conditions ;
496  val_b = 0 ; pari = 1 ;
497  for (int i=0; i<nr_lim; i++) {
498  val_b += pari*coef(lz, k, j, i) ;
499  pari = -pari ;
500  }
501  rhs.set(index) += val_b ;
502  index++ ;
503  if (derivative) {
504  der_b = 0 ; pari = -1 ;
505  for (int i=0; i<nr_lim; i++) {
506  der_b += pari*i*i*coef(lz, k, j, i) ;
507  pari = -pari ;
508  }
509  rhs.set(index) += der_b/alpha ;
510  index++ ;
511  }
512  val_b = 0 ;
513  for (int i=0; i<nr_lim; i++)
514  val_b += coef(lz, k, j, i) ;
515  der_b = 0 ;
516  for (int i=0; i<nr_lim; i++) {
517  der_b += i*i*coef(lz, k, j, i) ;
518  }
519  if (domain_bc==lz) {
520  rhs.set(index) = mu(k,j) - c_val*val_b - c_der*der_b/alpha ;
521  index++ ;
522  }
523  else {
524  rhs.set(index) = -val_b ;
525  if (derivative) rhs.set(index+1) = -der_b / alpha ;
526  }
527  }
528  Tbl solut(sys_size);
529  assert(l_q>=0) ;
530  switch (l_q) {
531  case 0 :
532  solut = sys_l0.inverse(rhs) ;
533  break ;
534  case 1:
535  solut = sys_l1.inverse(rhs) ;
536  break ;
537  default:
538  solut = (l_q%2 == 0 ? sys_even.inverse(rhs) :
539  sys_odd.inverse(rhs)) ;
540  }
541  index = 0 ;
542  int dec = (base_r == R_CHEBP ? 0 : 1) ;
543  nr = mgrid.get_nr(0) ;
544  if (isexcised==false){
545  if (l_q>=2) {
546  coef.set(0, k, j, nr-2-dec) = solut(index) ;
547  index++ ;
548  }
549  coef.set(0, k, j, nr-1-dec) = solut(index) ;
550  index++ ;
551  if (base_r == R_CHEBI)
552  coef.set(0, k, j, nr-1) = 0 ;
553  if (domain_bc > 0) {
554  //Pour domaine 1
555  nr = mgrid.get_nr(1) ;
556  for (nl=1; nl<=n_conditions; nl++) {
557  int ii = n_conditions - nl + 1 ;
558  coef.set(1, k, j, nr-ii) = solut(index) ;
559  index++ ;
560  }
561  }
562  }
563  else{
564  coef.set(1,k,j,nr-1)=solut(index);
565  index++;
566  }
567  for (int lz=2; lz<=domain_bc; lz++) {
568  nr = mgrid.get_nr(lz) ;
569  for (nl=1; nl<=n_conditions; nl++) {
570  int ii = n_conditions - nl + 1 ;
571  coef.set(lz, k, j, nr-ii) = solut(index) ;
572  index++ ;
573  }
574  }
575  } //End of nullite_plm
576  else {
577  for (int lz=0; lz<=domain_bc; lz++)
578  for (int i=0; i<mgrid.get_nr(lz); i++)
579  coef.set(lz, k, j, i) = 0 ;
580  }
581  } //End of loop on j
582  } //End of loop on k
583 }
584 
585 }
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::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::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:281
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:69
Lorene::Scalar::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:386
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::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::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:161
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::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