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
27char 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
65namespace Lorene {
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++) {
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)
134
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
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) {
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) ) ;
353 par_mat->get_matrice_mod(2) ) ;
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 }
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}
Bases of the spectral expansions.
Definition base_val.h:322
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Affine radial mapping.
Definition map.h:2027
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition map_af.C:477
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Matrix handling.
Definition matrice.h:152
Multi-domain grid.
Definition grilles.h:273
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:294
Parameter storage.
Definition param.h:125
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary).
Definition scalar.h:396
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va
Definition scalar.h:1294
void match_tau(Param &par_bc, Param *par_mat=0x0)
Method for matching accross domains and imposing boundary condition.
Valeur va
The numerical value of the Scalar
Definition scalar.h:405
Basic array class.
Definition tbl.h:161
int get_etat() const
Gives the logical state.
Definition tbl.h:394
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:372
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Mtbl * c
Values of the function at the points of the multi-grid
Definition valeur.h:299
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
void coef() const
Computes the coeffcients of *this.
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define R_CHEBP
base de Cheb. paire (rare) seulement
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Lorene prototypes.
Definition app_hor.h:64