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
68namespace Lorene {
69void 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::Base_val
Bases of the spectral expansions.
Definition base_val.h:325
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::Map_af
Affine radial mapping.
Definition map.h:2042
Lorene::Map_af::get_alpha
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition map_af.C:604
Lorene::Matrice
Matrix handling.
Definition matrice.h:152
Lorene::Matrice::set
double & set(int j, int i)
Read/write of a particuliar element.
Definition matrice.h:277
Lorene::Matrice::inverse
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition matrice.C:427
Lorene::Matrice::annule_hard
void annule_hard()
Sets the logical state to ETATQCQ (undefined state).
Definition matrice.C:196
Lorene::Matrice::set_lu
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition matrice.C:395
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:279
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::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_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:465
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::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::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:196
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::Param
Parameter storage.
Definition param.h:125
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:862
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
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:587
Lorene::Param::get_n_double_mod
int get_n_double_mod() const
Returns the number of stored modifiable double 's addresses.
Definition param.C:449
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::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::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::Param::get_n_int
int get_n_int() const
Returns the number of stored int 's addresses.
Definition param.C:242
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::Scalar::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:386
Lorene::Scalar::etat
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary).
Definition scalar.h:402
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::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::va
Valeur va
The numerical value of the Scalar.
Definition scalar.h:411
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::get_etat
int get_etat() const
Gives the logical state.
Definition tbl.h:394
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:375
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
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::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:301
Lorene
Lorene prototypes.
Definition app_hor.h:67