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