LORENE
sym_tensor_trans_dirac_bound2.C
1/*
2 * Methods to impose the Dirac gauge: divergence-free condition, with interior boundary conditions added
3 *
4 * (see file sym_tensor.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2006, 2007 Jerome Novak,Nicolas Vasset
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28// char sym_tensor_trans_dirac_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_trans_dirac_bound2.C,v 1.5 2015/08/10 15:32:26 j_novak Exp $" ;
29
30/*
31 * $Id: sym_tensor_trans_dirac_bound2.C,v 1.5 2015/08/10 15:32:26 j_novak Exp $
32 * $Log: sym_tensor_trans_dirac_bound2.C,v $
33 * Revision 1.5 2015/08/10 15:32:26 j_novak
34 * Better calls to Param::add_int(), to avoid weird problems (e.g. with g++ 4.8).
35 *
36 * Revision 1.4 2014/10/13 08:53:43 j_novak
37 * Lorene classes and functions now belong to the namespace Lorene.
38 *
39 * Revision 1.3 2014/10/06 15:13:19 j_novak
40 * Modified #include directives to use c++ syntax.
41 *
42 * Revision 1.2 2008/08/20 15:08:06 n_vasset
43 * Cleaning up the code...
44 *
45 * Revision 1.1 2007/05/04 16:45:25 n_vasset
46 * First version
47 *
48 * Revision 1.2 2006/10/24 13:03:19 j_novak
49 * New methods for the solution of the tensor wave equation. Perhaps, first
50 * operational version...
51 *
52 * Revision 1.1 2006/09/05 15:38:45 j_novak
53 * The fuctions sol_Dirac... are in a separate file, with new parameters to
54 * control the boundary conditions.
55 *
56 *
57 * $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_trans_dirac_bound2.C,v 1.5 2015/08/10 15:32:26 j_novak Exp $
58 *
59 */
60
61// C headers
62#include <cassert>
63#include <cmath>
64
65// Lorene headers
66#include "param.h"
67#include "param_elliptic.h"
68#include "metric.h"
69#include "diff.h"
70#include "proto.h"
71#include "param.h"
72
73
74//----------------------------------------------------------------------------------
75//
76// sol_Dirac_A
77//
78//----------------------------------------------------------------------------------
79
80namespace Lorene {
82
83 const Map_af* mp_aff = dynamic_cast<const Map_af*>(mp) ;
84 assert(mp_aff != 0x0) ; //Only affine mapping for the moment
85
86 // WARNING!
87 // This will only work if we have at least 2 shells!
88
89
90 const Mg3d& mgrid = *mp_aff->get_mg() ;
91 assert(mgrid.get_type_r(0) == RARE) ;
92 if (aaa.get_etat() == ETATZERO) {
93 tilde_mu = 0 ;
94 x_new = 0 ;
95 return ;
96 }
97
98 double dir = 1.;
99 double neum = 0.;
100
101 // On suppose que le sym_tensor_entré est bien transverse;
102
103 // assert ( maxabs(contract((*this).derive_cov(mets), 1, 2)) < 0.00000000000001 ) ;
104
105
106 bound_mu.set_spectral_va().ylm();
107
108
109 assert(aaa.get_etat() != ETATNONDEF) ;
110 int nz = mgrid.get_nzone() ;
111 int nzm1 = nz - 1 ;
112 bool ced = (mgrid.get_type_r(nzm1) == UNSURR) ;
113 int n_shell = ced ? nz-2 : nzm1 ;
114 int nz_bc = nzm1 ;
115 if (par_bc != 0x0)
116 if (par_bc->get_n_int() > 0) nz_bc = par_bc->get_int() ;
117 n_shell = (nz_bc < n_shell ? nz_bc : n_shell) ;
118 bool cedbc = (ced && (nz_bc == nzm1)) ;
119#ifndef NDEBUG
120 if (!cedbc) {
121 assert(par_bc != 0x0) ;
122 assert(par_bc->get_n_tbl_mod() >= 3) ;
123 }
124#endif
125 int nt = mgrid.get_nt(0) ;
126 int np = mgrid.get_np(0) ;
127
128 Scalar source = aaa ;
130 source_coq.annule_domain(0) ;
131 if (ced) source_coq.annule_domain(nzm1) ;
132 source_coq.mult_r() ;
133 source.set_spectral_va().ylm() ;
134 source_coq.set_spectral_va().ylm() ;
135 Base_val base = source.get_spectral_base() ;
136 base.mult_x() ;
137
138 tilde_mu.annule_hard() ;
139 tilde_mu.set_spectral_base(base) ;
140 tilde_mu.set_spectral_va().set_etat_cf_qcq() ;
141 tilde_mu.set_spectral_va().c_cf->annule_hard() ;
142 x_new.annule_hard() ;
143 x_new.set_spectral_base(base) ;
144 x_new.set_spectral_va().set_etat_cf_qcq() ;
145 x_new.set_spectral_va().c_cf->annule_hard() ;
146
147 Mtbl_cf sol_part_mu(mgrid, base) ; sol_part_mu.annule_hard() ;
148 Mtbl_cf sol_part_x(mgrid, base) ; sol_part_x.annule_hard() ;
149 Mtbl_cf sol_hom1_mu(mgrid, base) ; sol_hom1_mu.annule_hard() ;
150 Mtbl_cf sol_hom1_x(mgrid, base) ; sol_hom1_x.annule_hard() ;
151 Mtbl_cf sol_hom2_mu(mgrid, base) ; sol_hom2_mu.annule_hard() ;
152 Mtbl_cf sol_hom2_x(mgrid, base) ; sol_hom2_x.annule_hard() ;
153
154 int l_q, m_q, base_r ;
155
156 //---------------
157 //-- NUCLEUS ---
158 //---------------
159
160 // On va annuler toutes les solutions dans le noyau
161
162 { int lz = 0 ;
163 int nr = mgrid.get_nr(lz) ;
164 // double alpha = mp_aff->get_alpha()[lz] ;
165 // Matrice ope(2*nr, 2*nr) ;
166 // ope.set_etat_qcq() ;
167
168 for (int k=0 ; k<np+1 ; k++) {
169 for (int j=0 ; j<nt ; j++) {
170 // quantic numbers and spectral bases
171 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
172 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
173
174
175 Tbl sol(2*nr) ;
176 sol.set_etat_zero();
177 for (int i=0; i<nr; i++) {
178 sol_part_mu.set(lz, k, j, i) = 0 ;
179 sol_part_x.set(lz, k, j, i) = 0 ;
180 }
181 for (int i=0; i<nr; i++) {
182 sol_hom2_mu.set(lz, k, j, i) = 0 ;
183 sol_hom2_x.set(lz, k, j, i) = 0 ;
184 }
185 }
186 }
187 }
188 }
189
190 //-------------
191 // -- Shells --
192 //-------------
193
194 for (int lz=1; lz <= n_shell; lz++) {
195 int nr = mgrid.get_nr(lz) ;
196 assert(mgrid.get_nt(lz) == nt) ;
197 assert(mgrid.get_np(lz) == np) ;
198 double alpha = mp_aff->get_alpha()[lz] ;
199 double ech = mp_aff->get_beta()[lz] / alpha ;
200 Matrice ope(2*nr, 2*nr) ;
201 ope.set_etat_qcq() ;
202
203 for (int k=0 ; k<np+1 ; k++) {
204 for (int j=0 ; j<nt ; j++) {
205 // quantic numbers and spectral bases
206 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
207 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
208 Diff_xdsdx oxd(base_r, nr) ; const Matrice& mxd = oxd.get_matrice() ;
209 Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
210 Diff_id oid(base_r, nr) ; const Matrice& mid = oid.get_matrice() ;
211
212 for (int lin=0; lin<nr; lin++)
213 for (int col=0; col<nr; col++)
214 ope.set(lin,col) = mxd(lin,col) + ech*md(lin,col)
215 + 3*mid(lin,col) ;
216 for (int lin=0; lin<nr; lin++)
217 for (int col=0; col<nr; col++)
218 ope.set(lin,col+nr) = (2-l_q*(l_q+1))*mid(lin,col) ;
219 for (int lin=0; lin<nr; lin++)
220 for (int col=0; col<nr; col++)
221 ope.set(lin+nr,col) = -mid(lin,col) ;
222 for (int lin=0; lin<nr; lin++)
223 for (int col=0; col<nr; col++)
224 ope.set(lin+nr,col+nr) = mxd(lin,col) + ech*md(lin,col) ;
225
226 int ind0 = 0 ;
227 int ind1 = nr ;
228 for (int col=0; col<2*nr; col++) {
229 ope.set(ind0+nr-1, col) = 0 ;
230 ope.set(ind1+nr-1, col) = 0 ;
231 }
232 ope.set(ind0+nr-1, ind0) = 1 ;
233 ope.set(ind1+nr-1, ind1) = 1 ;
234
235 ope.set_lu() ;
236
237 Tbl sec(2*nr) ;
238 sec.set_etat_qcq() ;
239 for (int lin=0; lin<nr; lin++)
240 sec.set(lin) = 0 ;
241 for (int lin=0; lin<nr; lin++)
242 sec.set(nr+lin) = (*source_coq.get_spectral_va().c_cf)
243 (lz, k, j, lin) ;
244 sec.set(ind0+nr-1) = 0 ;
245 sec.set(ind1+nr-1) = 0 ;
246 Tbl sol = ope.inverse(sec) ;
247 for (int i=0; i<nr; i++) {
248 sol_part_mu.set(lz, k, j, i) = sol(i) ;
249 sol_part_x.set(lz, k, j, i) = sol(i+nr) ;
250 }
251 sec.annule_hard() ;
252 sec.set(ind0+nr-1) = 1 ;
253 sol = ope.inverse(sec) ;
254 for (int i=0; i<nr; i++) {
255 sol_hom1_mu.set(lz, k, j, i) = sol(i) ;
256 sol_hom1_x.set(lz, k, j, i) = sol(i+nr) ;
257 }
258 sec.set(ind0+nr-1) = 0 ;
259 sec.set(ind1+nr-1) = 1 ;
260 sol = ope.inverse(sec) ;
261 for (int i=0; i<nr; i++) {
262 sol_hom2_mu.set(lz, k, j, i) = sol(i) ;
263 sol_hom2_x.set(lz, k, j, i) = sol(i+nr) ;
264 }
265 }
266 }
267 }
268 }
269
270 //------------------------------
271 // Compactified external domain
272 //------------------------------
273 if (cedbc) {int lz = nzm1 ;
274 int nr = mgrid.get_nr(lz) ;
275 assert(mgrid.get_nt(lz) == nt) ;
276 assert(mgrid.get_np(lz) == np) ;
277 double alpha = mp_aff->get_alpha()[lz] ;
278 Matrice ope(2*nr, 2*nr) ;
279 ope.set_etat_qcq() ;
280
281 for (int k=0 ; k<np+1 ; k++) {
282 for (int j=0 ; j<nt ; j++) {
283 // quantic numbers and spectral bases
284 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
285 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
286 Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
287 Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
288
289 for (int lin=0; lin<nr; lin++)
290 for (int col=0; col<nr; col++)
291 ope.set(lin,col) = - md(lin,col) + 3*ms(lin,col) ;
292 for (int lin=0; lin<nr; lin++)
293 for (int col=0; col<nr; col++)
294 ope.set(lin,col+nr) = (2-l_q*(l_q+1))*ms(lin,col) ;
295 for (int lin=0; lin<nr; lin++)
296 for (int col=0; col<nr; col++)
297 ope.set(lin+nr,col) = -ms(lin,col) ;
298 for (int lin=0; lin<nr; lin++)
299 for (int col=0; col<nr; col++)
300 ope.set(lin+nr,col+nr) = -md(lin,col) ;
301
302 ope *= 1./alpha ;
303 int ind0 = 0 ;
304 int ind1 = nr ;
305 for (int col=0; col<2*nr; col++) {
306 ope.set(ind0+nr-1, col) = 0 ;
307 ope.set(ind1+nr-2, col) = 0 ;
308 ope.set(ind1+nr-1, col) = 0 ;
309 }
310 for (int col=0; col<nr; col++) {
311 ope.set(ind0+nr-1, col+ind0) = 1 ;
312 ope.set(ind1+nr-1, col+ind1) = 1 ;
313 }
314 ope.set(ind1+nr-2, ind1+1) = 1 ;
315
316 ope.set_lu() ;
317
318 Tbl sec(2*nr) ;
319 sec.set_etat_qcq() ;
320 for (int lin=0; lin<nr; lin++)
321 sec.set(lin) = 0 ;
322 for (int lin=0; lin<nr; lin++)
323 sec.set(nr+lin) = (*source.get_spectral_va().c_cf)
324 (lz, k, j, lin) ;
325 sec.set(ind0+nr-1) = 0 ;
326 sec.set(ind1+nr-2) = 0 ;
327 sec.set(ind1+nr-1) = 0 ;
328 Tbl sol = ope.inverse(sec) ;
329 for (int i=0; i<nr; i++) {
330 sol_part_mu.set(lz, k, j, i) = sol(i) ;
331 sol_part_x.set(lz, k, j, i) = sol(i+nr) ;
332 }
333 sec.annule_hard() ;
334 sec.set(ind1+nr-2) = 1 ;
335 sol = ope.inverse(sec) ;
336 for (int i=0; i<nr; i++) {
337 sol_hom1_mu.set(lz, k, j, i) = sol(i) ;
338 sol_hom1_x.set(lz, k, j, i) = sol(i+nr) ;
339 }
340 }
341 }
342 }
343 }
344
345
346 // On écrit maintenant le systčme de raccord
347
348 int taille = 2*nz_bc ;
349 if (cedbc) taille-- ;
350 Mtbl_cf& mmu = *tilde_mu.set_spectral_va().c_cf ;
351 Mtbl_cf& mw = *x_new.set_spectral_va().c_cf ;
352
353 Tbl sec_membre(taille) ;
354 Matrice systeme(taille, taille) ;
355 int ligne ; int colonne ;
356 Tbl pipo(1) ;
357 const Tbl& mub = (cedbc ? pipo : par_bc->get_tbl_mod(2) );
358 double c_mu = (cedbc ? 0 : par_bc->get_tbl_mod(0)(0) ) ;
359 double d_mu = (cedbc ? 0 : par_bc->get_tbl_mod(0)(1) ) ;
360 double c_x = (cedbc ? 0 : par_bc->get_tbl_mod(0)(2) ) ;
361 double d_x = (cedbc ? 0 : par_bc->get_tbl_mod(0)(3) ) ;
368
369
370 dhom1_mu.dsdx() ;
371 dhom2_mu.dsdx() ;
372 dpart_mu.dsdx() ;
373 dhom1_x.dsdx() ;
374 dhom2_x.dsdx() ;
375 dpart_x.dsdx() ;
376
377
378 // Loop on l and m
379 //----------------
380 for (int k=0 ; k<np+1 ; k++)
381 for (int j=0 ; j<nt ; j++) {
382 base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
383 if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
384 ligne = 0 ;
385 colonne = 0 ;
386 systeme.annule_hard() ;
387 sec_membre.annule_hard() ;
388
389 //First shell
390 // Internal boundary condition
391 int nr = mgrid.get_nr(1) ;
392
393 systeme.set(ligne, colonne) = dir * sol_hom1_mu.val_in_bound_jk(1, j, k) + neum* dhom1_mu.val_in_bound_jk(1,j,k);
394 systeme.set(ligne, colonne+1) = dir * sol_hom2_mu.val_in_bound_jk(1, j, k) + neum* dhom2_mu.val_in_bound_jk(1,j,k);
395
396 Mtbl_cf *boundmu = bound_mu.get_spectral_va().c_cf;
397 double blob = (*boundmu).val_in_bound_jk(1,j,k);
398 sec_membre.set(ligne) = - dir* sol_part_mu.val_in_bound_jk(1, j, k) - neum * dpart_mu.val_in_bound_jk(1,j,k) + blob;
399
400
401
402 ligne++;
403 // Condition at x=1
404 systeme.set(ligne, colonne) =
405 sol_hom1_mu.val_out_bound_jk(1, j, k) ;
406 systeme.set(ligne, colonne+1) =
407 sol_hom2_mu.val_out_bound_jk(1, j, k) ;
408
409 sec_membre.set(ligne) -= sol_part_mu.val_out_bound_jk(1, j, k) ;
410 ligne++ ;
411
412 systeme.set(ligne, colonne) =
413 sol_hom1_x.val_out_bound_jk(1, j, k) ;
414 systeme.set(ligne, colonne+1) =
415 sol_hom2_x.val_out_bound_jk(1, j, k) ;
416
417 sec_membre.set(ligne) -= sol_part_x.val_out_bound_jk(1, j, k) ;
418
419 colonne += 2 ;
420
421 //shells with nz>2
422 for (int zone=2 ; zone<nz_bc ; zone++) {
423 nr = mgrid.get_nr(zone) ;
424 ligne-- ;
425
426 //Condition at x = -1
427 systeme.set(ligne, colonne) =
428 - sol_hom1_mu.val_in_bound_jk(zone, j, k) ;
429 systeme.set(ligne, colonne+1) =
430 - sol_hom2_mu.val_in_bound_jk(zone, j, k) ;
431
432 sec_membre.set(ligne) += sol_part_mu.val_in_bound_jk(zone, j, k) ;
433 ligne++ ;
434
435 systeme.set(ligne, colonne) =
436 - sol_hom1_x.val_in_bound_jk(zone, j, k) ;
437 systeme.set(ligne, colonne+1) =
438 - sol_hom2_x.val_in_bound_jk(zone, j, k) ;
439
440 sec_membre.set(ligne) += sol_part_x.val_in_bound_jk(zone, j, k) ;
441 ligne++ ;
442
443 // Condition at x=1
444 systeme.set(ligne, colonne) =
445 sol_hom1_mu.val_out_bound_jk(zone, j, k) ;
446 systeme.set(ligne, colonne+1) =
447 sol_hom2_mu.val_out_bound_jk(zone, j, k) ;
448
449 sec_membre.set(ligne) -= sol_part_mu.val_out_bound_jk(zone, j, k) ;
450 ligne++ ;
451
452 systeme.set(ligne, colonne) =
453 sol_hom1_x.val_out_bound_jk(zone, j, k) ;
454 systeme.set(ligne, colonne+1) =
455 sol_hom2_x.val_out_bound_jk(zone, j, k) ;
456
457 sec_membre.set(ligne) -= sol_part_x.val_out_bound_jk(zone, j, k) ;
458
459 colonne += 2 ;
460 }
461
462 //Last domain
463 nr = mgrid.get_nr(nz_bc) ;
464 double alpha = mp_aff->get_alpha()[nz_bc] ;
465 ligne-- ;
466
467 //Condition at x = -1
468 systeme.set(ligne, colonne) =
469 - sol_hom1_mu.val_in_bound_jk(nz_bc, j, k) ;
470 if (!cedbc) systeme.set(ligne, colonne+1) =
471 - sol_hom2_mu.val_in_bound_jk(nz_bc, j, k) ;
472
473 sec_membre.set(ligne) += sol_part_mu.val_in_bound_jk(nz_bc, j, k) ;
474 ligne++ ;
475
476 systeme.set(ligne, colonne) =
477 - sol_hom1_x.val_in_bound_jk(nz_bc, j, k) ;
478 if (!cedbc) systeme.set(ligne, colonne+1) =
479 - sol_hom2_x.val_in_bound_jk(nz_bc, j, k) ;
480
481 sec_membre.set(ligne) += sol_part_x.val_in_bound_jk(nz_bc, j, k) ;
482 ligne++ ;
483
484 if (!cedbc) {// Special condition at x=1
485
486 systeme.set(ligne, colonne) =
487 c_mu*sol_hom1_mu.val_out_bound_jk(nz_bc, j, k)
488 + d_mu*dhom1_mu.val_out_bound_jk(nz_bc, j, k) / alpha
489 + c_x*sol_hom1_x.val_out_bound_jk(nz_bc, j, k)
490 + d_x*dhom1_x.val_out_bound_jk(nz_bc, j, k) / alpha ;
491 systeme.set(ligne, colonne+1) =
492 c_mu*sol_hom2_mu.val_out_bound_jk(nz_bc, j, k)
493 + d_mu*dhom2_mu.val_out_bound_jk(nz_bc, j, k) / alpha
494 + c_x*sol_hom2_x.val_out_bound_jk(nz_bc, j, k)
495 + d_x*dhom2_x.val_out_bound_jk(nz_bc, j, k) / alpha ;
496
497 sec_membre.set(ligne) -= c_mu*sol_part_mu.val_out_bound_jk(nz_bc, j, k)
498 + d_mu*dpart_mu.val_out_bound_jk(nz_bc, j, k)/alpha
499 + c_x*sol_part_mu.val_out_bound_jk(nz_bc, j, k)
500 + d_x*dpart_mu.val_out_bound_jk(nz_bc, j, k)/alpha
501 - mub(k, j) ;
502 }
503
504
505
506
508 // System inversion: Solution of the system giving the coefficients for the homogeneous
509 // solutions
510 //-------------------------------------------------------------------
511 systeme.set_lu() ;
512 Tbl facteur = systeme.inverse(sec_membre) ;
513
514
515 int conte = 0 ;
516
517 // everything is put to the right place...
518 //----------------------------------------
519 nr = mgrid.get_nr(0) ; //nucleus
520 for (int i=0 ; i<nr ; i++) {
521 mmu.set(0, k, j, i) = 0 ;
522 mw.set(0, k, j, i) = 0 ;
523 }
524 nr = mgrid.get_nr(1) ; //first shell
525 for (int i=0 ; i<nr ; i++) {
526 mmu.set(1, k, j, i) = sol_part_mu(1, k, j, i)
527 + facteur(conte)*sol_hom1_mu(1, k, j, i)
528 + facteur(conte+1)*sol_hom2_mu(1, k, j, i);
529 mw.set(1, k, j, i) = sol_part_x(1, k, j, i)
530 + facteur(conte)*sol_hom1_x(1, k, j, i)
531 + facteur(conte+1)*sol_hom2_x(1,k,j,i);
532 }
533 conte+=2 ;
534 for (int zone=2 ; zone<=n_shell ; zone++) { //shells
535 nr = mgrid.get_nr(zone) ;
536 for (int i=0 ; i<nr ; i++) {
537 mmu.set(zone, k, j, i) = sol_part_mu(zone, k, j, i)
539 + facteur(conte+1)*sol_hom2_mu(zone, k, j, i) ;
540
541 mw.set(zone, k, j, i) = sol_part_x(zone, k, j, i)
543 + facteur(conte+1)*sol_hom2_x(zone, k, j, i) ;
544 }
545 conte+=2 ;
546 }
547 if (cedbc) {
548 nr = mgrid.get_nr(nzm1) ; //compactified external domain
549 for (int i=0 ; i<nr ; i++) {
550 mmu.set(nzm1, k, j, i) = sol_part_mu(nzm1, k, j, i)
551 + facteur(conte)*sol_hom1_mu(nzm1, k, j, i) ;
552
553 mw.set(nzm1, k, j, i) = sol_part_x(nzm1, k, j, i)
554 + facteur(conte)*sol_hom1_x(nzm1, k, j, i) ;
555 }
556 }
557 } // End of nullite_plm
558 } //End of loop on theta
559
560 if (tilde_mu.set_spectral_va().c != 0x0)
561 delete tilde_mu.set_spectral_va().c ;
562 tilde_mu.set_spectral_va().c = 0x0 ;
563 tilde_mu.set_spectral_va().ylm_i() ;
564
565 if (x_new.set_spectral_va().c != 0x0)
566 delete x_new.set_spectral_va().c ;
567 x_new.set_spectral_va().c = 0x0 ;
568 x_new.set_spectral_va().ylm_i();
569
570}
571
572
573
574//----------------------------------------------------------------------------------
575//
576// sol_Dirac_tilde_B
577//
578//----------------------------------------------------------------------------------
579
580void Sym_tensor_trans::sol_Dirac_BC2(const Scalar& bb, const Scalar& cc, const Scalar& hh,
581 Scalar& hrr, Scalar& tilde_eta, Scalar& ww, Scalar bound_eta,double dir, double neum, double rhor, Param* par_bc, Param* par_mat) {
582
583 const Map_af* mp_aff = dynamic_cast<const Map_af*>(mp) ;
584 assert(mp_aff != 0x0) ; //Only affine mapping for the moment
585
586 const Mg3d& mgrid = *mp_aff->get_mg() ;
587 assert(mgrid.get_type_r(0) == RARE) ;
588
589 int nz = mgrid.get_nzone() ;
590 int nzm1 = nz - 1 ;
591 bool ced = (mgrid.get_type_r(nzm1) == UNSURR) ;
592 int n_shell = ced ? nz-2 : nzm1 ;
593 int nz_bc = nzm1 ;
594 if (par_bc != 0x0)
595 if (par_bc->get_n_int() > 0)
596 nz_bc = par_bc->get_int() ;
597 n_shell = (nz_bc < n_shell ? nz_bc : n_shell) ;
598 bool cedbc = (ced && (nz_bc == nzm1)) ;
599#ifndef NDEBUG
600 if (!cedbc) {
601 assert(par_bc != 0x0) ;
602 assert(par_bc->get_n_tbl_mod() >= 2) ;
603 }
604#endif
605 int nt = mgrid.get_nt(0) ;
606 int np = mgrid.get_np(0) ;
607
608 int l_q, m_q, base_r;
609
610 // ON passe en ylm les quantités B et C !!! CRUCIAL !!!
611
612 Scalar bb2 = bb; bb2.set_spectral_va().ylm();
613 Scalar cc2 = cc; cc2.set_spectral_va().ylm();
614
615
616 Scalar hh2 = hh; hh2.set_spectral_va().ylm();
617
618 Base_val base = hh2.get_spectral_base() ;
619
620
621 Scalar tilde_b(*mp_aff); tilde_b.annule_hard(); tilde_b.std_spectral_base();
622 tilde_b.set_spectral_va().ylm();
623
624 // Base_val base = hrr.get_spectral_base();
625
626 // int m_q, l_q, base_r ;
627 for (int lz=0; lz<nz; lz++) {
628 int nr = mgrid.get_nr(lz);
629 for (int k=0; k<np+1; k++)
630 for (int j=0; j<nt; j++) {
631 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
632 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1))
633 {
634 for (int i=0; i<nr; i++)
635 tilde_b.set_spectral_va().c_cf->set(lz, k, j, i)
636 += 2*(*bb2.get_spectral_va().c_cf)(lz, k, j, i)
637 + (*cc2.get_spectral_va().c_cf)(lz, k, j, i)/ (2* double(l_q +1.));
638 }
639 }
640 }
641
642
643 if (tilde_b.set_spectral_va().c != 0x0)
644 delete tilde_b.set_spectral_va().c ;
645 tilde_b.set_spectral_va().c = 0x0 ;
646 tilde_b.set_spectral_va().coef_i();
647
648
649 if ( (tilde_b.get_etat() == ETATZERO) && (hh.get_etat() == ETATZERO) ) {
650 hrr = 0 ;
651 tilde_eta = 0 ;
652 ww = 0 ;
653 return ;
654 }
655
656 bound_eta.set_spectral_va().ylm();
658 tilde_b2.set_spectral_va().ylm();
659
660 assert (tilde_b.get_etat() != ETATNONDEF) ;
661 assert (hh.get_etat() != ETATNONDEF) ;
662
665 source_coq.annule_domain(0) ;
666 source_coq.annule_domain(nzm1) ;
667 source_coq.mult_r() ;
668 source.set_spectral_va().ylm() ;
669 source_coq.set_spectral_va().ylm() ;
670 bool bnull = (tilde_b.get_etat() == ETATZERO) ;
671
672 assert(hh.check_dzpuis(0)) ;
673 Scalar hoverr = hh ;
674 hoverr.div_r_dzpuis(2) ;
675 hoverr.set_spectral_va().ylm() ;
676 Scalar dhdr = hh.dsdr() ;
677 dhdr.set_spectral_va().ylm() ;
678 Scalar h_coq = hh ;
679 h_coq.set_spectral_va().ylm() ;
680 Scalar dh_coq = hh.dsdr() ;
681 dh_coq.mult_r_dzpuis(0) ;
682 dh_coq.set_spectral_va().ylm() ;
683 bool hnull = (hh.get_etat() == ETATZERO) ;
684
685 int lmax = base.give_lmax(mgrid, 0) + 1;
686
687
688 // Utilisation des params, facultative, permettant uniquement une optimisation....
689
690 bool need_calculation = true ;
691 if (par_mat != 0x0) {
692 bool param_new = false ;
693 if ((par_mat->get_n_int_mod() >= 4)
694 &&(par_mat->get_n_tbl_mod()>=1)
695 &&(par_mat->get_n_matrice_mod()>=1)
696 &&(par_mat->get_n_itbl_mod()>=1)) {
697 if (par_mat->get_int_mod(0) < nz_bc) param_new = true ;
698 if (par_mat->get_int_mod(1) != lmax) param_new = true ;
699 if (par_mat->get_int_mod(2) != mgrid.get_type_t() ) param_new = true ;
700 if (par_mat->get_int_mod(3) != mgrid.get_type_p() ) param_new = true ;
701 if (par_mat->get_itbl_mod(0)(0) != mgrid.get_nr(0)) param_new = true ;
702 if (fabs(par_mat->get_tbl_mod(0)(0) - mp_aff->get_alpha()[0]) > 2.e-15)
703 param_new = true ;
704 for (int l=1; l<= n_shell; l++) {
705 if (par_mat->get_itbl_mod(0)(l) != mgrid.get_nr(l)) param_new = true ;
706 if (fabs(par_mat->get_tbl_mod(0)(l) - mp_aff->get_beta()[l] /
707 mp_aff->get_alpha()[l]) > 2.e-15) param_new = true ;
708 }
709 if (ced) {
710 if (par_mat->get_itbl_mod(0)(nzm1) != mgrid.get_nr(nzm1)) param_new = true ;
711 if (fabs(par_mat->get_tbl_mod(0)(nzm1) - mp_aff->get_alpha()[nzm1]) > 2.e-15)
712 param_new = true ;
713 }
714 }
715 else{
716 param_new = true ;
717 }
718 if (param_new) {
719 par_mat->clean_all() ;
720 int* nz_bc_new = new int(nz_bc) ;
721 par_mat->add_int_mod(*nz_bc_new, 0) ;
722 int* lmax_new = new int(lmax) ;
723 par_mat->add_int_mod(*lmax_new, 1) ;
724 int* type_t_new = new int(mgrid.get_type_t()) ;
725 par_mat->add_int_mod(*type_t_new, 2) ;
726 int* type_p_new = new int(mgrid.get_type_p()) ;
727 par_mat->add_int_mod(*type_p_new, 3) ;
728 Itbl* pnr = new Itbl(nz) ;
729 pnr->set_etat_qcq() ;
730 par_mat->add_itbl_mod(*pnr) ;
731 for (int l=0; l<nz; l++)
732 pnr->set(l) = mgrid.get_nr(l) ;
733 Tbl* palpha = new Tbl(nz) ;
734 palpha->set_etat_qcq() ;
735 par_mat->add_tbl_mod(*palpha) ;
736 palpha->set(0) = mp_aff->get_alpha()[0] ;
737 for (int l=1; l<nzm1; l++)
738 palpha->set(l) = mp_aff->get_beta()[l] / mp_aff->get_alpha()[l] ;
739 palpha->set(nzm1) = mp_aff->get_alpha()[nzm1] ;
740 }
741 else need_calculation = false ;
742 }
743
744 hrr.set_etat_qcq() ;
745 hrr.set_spectral_base(base) ;
746 hrr.set_spectral_va().set_etat_cf_qcq() ;
747 hrr.set_spectral_va().c_cf->annule_hard() ;
748 tilde_eta.annule_hard() ;
749 tilde_eta.set_spectral_base(base) ;
750 tilde_eta.set_spectral_va().set_etat_cf_qcq() ;
751 tilde_eta.set_spectral_va().c_cf->annule_hard() ;
752 ww.annule_hard() ;
753 ww.set_spectral_base(base) ;
756
758 tilde_eta.annule_l(0,0, true) ;
759
760 Mtbl_cf sol_part_hrr(mgrid, base) ; sol_part_hrr.annule_hard() ;
761 Mtbl_cf sol_part_eta(mgrid, base) ; sol_part_eta.annule_hard() ;
762 Mtbl_cf sol_part_w(mgrid, base) ; sol_part_w.annule_hard() ;
763 Mtbl_cf sol_hom1_hrr(mgrid, base) ; sol_hom1_hrr.annule_hard() ;
764 Mtbl_cf sol_hom1_eta(mgrid, base) ; sol_hom1_eta.annule_hard() ;
765 Mtbl_cf sol_hom1_w(mgrid, base) ; sol_hom1_w.annule_hard() ;
766 Mtbl_cf sol_hom2_hrr(mgrid, base) ; sol_hom2_hrr.annule_hard() ;
767 Mtbl_cf sol_hom2_eta(mgrid, base) ; sol_hom2_eta.annule_hard() ;
768 Mtbl_cf sol_hom2_w(mgrid, base) ; sol_hom2_w.annule_hard() ;
769 Mtbl_cf sol_hom3_hrr(mgrid, base) ; sol_hom3_hrr.annule_hard() ;
770 Mtbl_cf sol_hom3_eta(mgrid, base) ; sol_hom3_eta.annule_hard() ;
771 Mtbl_cf sol_hom3_w(mgrid, base) ; sol_hom3_w.annule_hard() ;
772
774
775
776 // Calcul des solutions homogenes et particulieres...
777
778
779 //---------------
780 //-- NUCLEUS ---
781 //---------------
782 {int lz = 0 ;
783 int nr = mgrid.get_nr(lz) ;
784 double alpha = mp_aff->get_alpha()[lz] ;
785 Matrice ope(3*nr, 3*nr) ;
786 int ind2 = 2*nr ;
787 if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
788
789 for (int k=0 ; k<np+1 ; k++) {
790 for (int j=0 ; j<nt ; j++) {
791 // quantic numbers and spectral bases
792 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
793 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
794 if (need_calculation) {
795 ope.set_etat_qcq() ;
796 Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
797 Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
798
799 for (int lin=0; lin<nr; lin++)
800 for (int col=0; col<nr; col++)
801 ope.set(lin,col) = md(lin,col) + 3*ms(lin,col) ;
802 for (int lin=0; lin<nr; lin++)
803 for (int col=0; col<nr; col++)
804 ope.set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
805 for (int lin=0; lin<nr; lin++)
806 for (int col=0; col<nr; col++)
807 ope.set(lin,col+2*nr) = 0 ;
808 for (int lin=0; lin<nr; lin++)
809 for (int col=0; col<nr; col++)
810 ope.set(lin+nr,col) = -0.5*ms(lin,col) ;
811 for (int lin=0; lin<nr; lin++)
812 for (int col=0; col<nr; col++)
813 ope.set(lin+nr,col+nr) = md(lin,col) + 3*ms(lin,col) ;
814 for (int lin=0; lin<nr; lin++)
815 for (int col=0; col<nr; col++)
816 ope.set(lin+nr,col+2*nr) = (2. - l_q*(l_q+1))*ms(lin,col) ;
817 for (int lin=0; lin<nr; lin++)
818 for (int col=0; col<nr; col++)
819 ope.set(lin+2*nr,col) = -0.5*md(lin,col)/double(l_q+1)
820 - 0.5*double(l_q+4)/double(l_q+1)*ms(lin,col) ;
821 for (int lin=0; lin<nr; lin++)
822 for (int col=0; col<nr; col++)
823 ope.set(lin+2*nr,col+nr) = -2*ms(lin,col) ;
824 for (int lin=0; lin<nr; lin++)
825 for (int col=0; col<nr; col++)
826 ope.set(lin+2*nr,col+2*nr) = (l_q+2)*md(lin,col)
827 + l_q*(l_q+2)*ms(lin,col) ;
828
829 ope *= 1./alpha ;
830 for (int col=0; col<3*nr; col++)
831 if (l_q>2) ope.set(ind2+nr-2, col) = 0 ;
832 for (int col=0; col<3*nr; col++) {
833 ope.set(nr-1, col) = 0 ;
834 ope.set(2*nr-1, col) = 0 ;
835 ope.set(3*nr-1, col) = 0 ;
836 }
837 int pari = 1 ;
838 if (base_r == R_CHEBP) {
839 for (int col=0; col<nr; col++) {
840 ope.set(nr-1, col) = pari ;
841 ope.set(2*nr-1, col+nr) = pari ;
842 ope.set(3*nr-1, col+2*nr) = pari ;
843 pari = - pari ;
844 }
845 }
846 else { //In the odd case, the last coefficient must be zero!
847 ope.set(nr-1, nr-1) = 1 ;
848 ope.set(2*nr-1, 2*nr-1) = 1 ;
849 ope.set(3*nr-1, 3*nr-1) = 1 ;
850 }
851 if (l_q>2)
852 ope.set(ind2+nr-2, ind2) = 1 ;
853
854 ope.set_lu() ;
855 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
856 Matrice* pope = new Matrice(ope) ;
857 par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
858 mat_done.set(l_q) = 1 ;
859 }
860 } //End of case when a calculation is needed
861
862 const Matrice& oper = (par_mat == 0x0 ? ope :
863 par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
864 Tbl sec(3*nr) ;
865 sec.set_etat_qcq() ;
866 if (hnull) {
867 for (int lin=0; lin<2*nr; lin++)
868 sec.set(lin) = 0 ;
869 for (int lin=0; lin<nr; lin++)
870 sec.set(2*nr+lin) = (*source.get_spectral_va().c_cf)
871 (lz, k, j, lin) ;
872 }
873 else {
874 for (int lin=0; lin<nr; lin++)
875 sec.set(lin) = (*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) ;
876 for (int lin=0; lin<nr; lin++)
877 sec.set(lin+nr) = -0.5*(*hoverr.get_spectral_va().c_cf)
878 (lz, k, j, lin) ;
879 if (bnull) {
880 for (int lin=0; lin<nr; lin++)
881 sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
882 (*dhdr.get_spectral_va().c_cf)(lz, k, j, lin)
883 + (l_q+2)*(*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) ) ;
884 }
885 else {
886 for (int lin=0; lin<nr; lin++)
887 sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
888 (*dhdr.get_spectral_va().c_cf)(lz, k, j, lin)
889 + (l_q+2)*(*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) )
890 + (*source.get_spectral_va().c_cf)(lz, k, j, lin) ;
891 }
892 }
893 if (l_q>2) sec.set(ind2+nr-2) = 0 ;
894 sec.set(3*nr-1) = 0 ;
895 Tbl sol = oper.inverse(sec) ;
896 for (int i=0; i<nr; i++) {
897 sol_part_hrr.set(lz, k, j, i) = sol(i) ;
898 sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
899 sol_part_w.set(lz, k, j, i) = sol(i+2*nr) ;
900 }
901 sec.annule_hard() ;
902 if (l_q>2) {
903 sec.set(ind2+nr-2) = 1 ;
904 sol = oper.inverse(sec) ;
905 }
906 else { //Homogeneous solution put in by hand in the case l=2
907 sol.annule_hard() ;
908 sol.set(0) = 4 ;
909 sol.set(nr) = 2 ;
910 sol.set(2*nr) = 1 ;
911 }
912 for (int i=0; i<nr; i++) {
913 sol_hom3_hrr.set(lz, k, j, i) = sol(i) ;
914 sol_hom3_eta.set(lz, k, j, i) = sol(i+nr) ;
915 sol_hom3_w.set(lz, k, j, i) = sol(i+2*nr) ;
916 }
917 }
918 }
919 }
920 }
921
922
923 //-------------
924 // -- Shells --
925 //-------------
926
927 for (int lz=1; lz<= n_shell; lz++) {
928 if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
929 int nr = mgrid.get_nr(lz) ;
930 int ind0 = 0 ;
931 int ind1 = nr ;
932 int ind2 = 2*nr ;
933 double alpha = mp_aff->get_alpha()[lz] ;
934 double ech = mp_aff->get_beta()[lz] / alpha ;
935 Matrice ope(3*nr, 3*nr) ;
936
937 for (int k=0 ; k<np+1 ; k++) {
938 for (int j=0 ; j<nt ; j++) {
939 // quantic numbers and spectral bases
940 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
941 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
942 if (need_calculation) {
943 ope.set_etat_qcq() ;
944 Diff_xdsdx oxd(base_r, nr) ; const Matrice& mxd = oxd.get_matrice() ;
945 Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
946 Diff_id oid(base_r, nr) ; const Matrice& mid = oid.get_matrice() ;
947
948 for (int lin=0; lin<nr; lin++)
949 for (int col=0; col<nr; col++)
950 ope.set(lin,col) = mxd(lin,col) + ech*md(lin,col)
951 + 3*mid(lin,col) ;
952 for (int lin=0; lin<nr; lin++)
953 for (int col=0; col<nr; col++)
954 ope.set(lin,col+nr) = -l_q*(l_q+1)*mid(lin,col) ;
955 for (int lin=0; lin<nr; lin++)
956 for (int col=0; col<nr; col++)
957 ope.set(lin,col+2*nr) = 0 ;
958 for (int lin=0; lin<nr; lin++)
959 for (int col=0; col<nr; col++)
960 ope.set(lin+nr,col) = -0.5*mid(lin,col) ;
961 for (int lin=0; lin<nr; lin++)
962 for (int col=0; col<nr; col++)
963 ope.set(lin+nr,col+nr) = mxd(lin,col) + ech*md(lin,col)
964 + 3*mid(lin,col) ;
965 for (int lin=0; lin<nr; lin++)
966 for (int col=0; col<nr; col++)
967 ope.set(lin+nr,col+2*nr) = (2. - l_q*(l_q+1))*mid(lin,col) ;
968 for (int lin=0; lin<nr; lin++)
969 for (int col=0; col<nr; col++)
970 ope.set(lin+2*nr,col) =
971 -0.5/double(l_q+1)*(mxd(lin,col) + ech*md(lin,col)
972 + double(l_q+4)*mid(lin,col)) ;
973 for (int lin=0; lin<nr; lin++)
974 for (int col=0; col<nr; col++)
975 ope.set(lin+2*nr,col+nr) = -2*mid(lin,col) ;
976 for (int lin=0; lin<nr; lin++)
977 for (int col=0; col<nr; col++)
978 ope.set(lin+2*nr,col+2*nr) =
979 double(l_q+2)*(mxd(lin,col) + ech*md(lin,col)
980 + l_q*mid(lin,col)) ;
981 for (int col=0; col<3*nr; col++) {
982 ope.set(ind0+nr-1, col) = 0 ;
983 ope.set(ind1+nr-1, col) = 0 ;
984 ope.set(ind2+nr-1, col) = 0 ;
985 }
986 ope.set(ind0+nr-1, ind0) = 1 ;
987 ope.set(ind1+nr-1, ind1) = 1 ;
988 ope.set(ind2+nr-1, ind2) = 1 ;
989
990 ope.set_lu() ;
991 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
992 Matrice* pope = new Matrice(ope) ;
993 par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
994 mat_done.set(l_q) = 1 ;
995 }
996 } //End of case when a calculation is needed
997 const Matrice& oper = (par_mat == 0x0 ? ope :
998 par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
999 Tbl sec(3*nr) ;
1000 sec.set_etat_qcq() ;
1001 if (hnull) {
1002 for (int lin=0; lin<2*nr; lin++)
1003 sec.set(lin) = 0 ;
1004 for (int lin=0; lin<nr; lin++)
1005 sec.set(2*nr+lin) = (*source_coq.get_spectral_va().c_cf)
1006 (lz, k, j, lin) ;
1007 }
1008 else {
1009 for (int lin=0; lin<nr; lin++)
1010 sec.set(lin) = (*h_coq.get_spectral_va().c_cf)(lz, k, j, lin) ;
1011 for (int lin=0; lin<nr; lin++)
1012 sec.set(lin+nr) = -0.5*(*h_coq.get_spectral_va().c_cf)
1013 (lz, k, j, lin) ;
1014 if (bnull) {
1015 for (int lin=0; lin<nr; lin++)
1016 sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
1017 (*dh_coq.get_spectral_va().c_cf)(lz, k, j, lin)
1018 + (l_q+2)*(*h_coq.get_spectral_va().c_cf)(lz, k, j, lin) ) ;
1019 }
1020 else {
1021 for (int lin=0; lin<nr; lin++)
1022 sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
1023 (*dh_coq.get_spectral_va().c_cf)(lz, k, j, lin)
1024 + (l_q+2)*(*h_coq.get_spectral_va().c_cf)(lz, k, j, lin) )
1025 + (*source_coq.get_spectral_va().c_cf)(lz, k, j, lin) ;
1026 }
1027 }
1028 sec.set(ind0+nr-1) = 0 ;
1029 sec.set(ind1+nr-1) = 0 ;
1030 sec.set(ind2+nr-1) = 0 ;
1031 Tbl sol = oper.inverse(sec) ;
1032 for (int i=0; i<nr; i++) {
1033 sol_part_hrr.set(lz, k, j, i) = sol(i) ;
1034 sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
1035 sol_part_w.set(lz, k, j, i) = sol(i+2*nr) ;
1036 }
1037 sec.annule_hard() ;
1038 sec.set(ind0+nr-1) = 1 ;
1039 sol = oper.inverse(sec) ;
1040 for (int i=0; i<nr; i++) {
1041 sol_hom1_hrr.set(lz, k, j, i) = sol(i) ;
1042 sol_hom1_eta.set(lz, k, j, i) = sol(i+nr) ;
1043 sol_hom1_w.set(lz, k, j, i) = sol(i+2*nr) ;
1044 }
1045 sec.set(ind0+nr-1) = 0 ;
1046 sec.set(ind1+nr-1) = 1 ;
1047 sol = oper.inverse(sec) ;
1048 for (int i=0; i<nr; i++) {
1049 sol_hom2_hrr.set(lz, k, j, i) = sol(i) ;
1050 sol_hom2_eta.set(lz, k, j, i) = sol(i+nr) ;
1051 sol_hom2_w.set(lz, k, j, i) = sol(i+2*nr) ;
1052 }
1053 sec.set(ind1+nr-1) = 0 ;
1054 sec.set(ind2+nr-1) = 1 ;
1055 sol = oper.inverse(sec) ;
1056 for (int i=0; i<nr; i++) {
1057 sol_hom3_hrr.set(lz, k, j, i) = sol(i) ;
1058 sol_hom3_eta.set(lz, k, j, i) = sol(i+nr) ;
1059 sol_hom3_w.set(lz, k, j, i) = sol(i+2*nr) ;
1060 }
1061 }
1062 }
1063 }
1064 }
1065
1066 //------------------------------
1067 // Compactified external domain
1068 //------------------------------
1069 if (cedbc) {int lz = nzm1 ;
1070 if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
1071 int nr = mgrid.get_nr(lz) ;
1072 int ind0 = 0 ;
1073 int ind1 = nr ;
1074 int ind2 = 2*nr ;
1075 double alpha = mp_aff->get_alpha()[lz] ;
1076 Matrice ope(3*nr, 3*nr) ;
1077
1078 for (int k=0 ; k<np+1 ; k++) {
1079 for (int j=0 ; j<nt ; j++) {
1080 // quantic numbers and spectral bases
1081 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
1082 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
1083 if (need_calculation) {
1084 ope.set_etat_qcq() ;
1085 Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
1086 Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
1087
1088 for (int lin=0; lin<nr; lin++)
1089 for (int col=0; col<nr; col++)
1090 ope.set(lin,col) = - md(lin,col) + 3*ms(lin,col) ;
1091 for (int lin=0; lin<nr; lin++)
1092 for (int col=0; col<nr; col++)
1093 ope.set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
1094 for (int lin=0; lin<nr; lin++)
1095 for (int col=0; col<nr; col++)
1096 ope.set(lin,col+2*nr) = 0 ;
1097 for (int lin=0; lin<nr; lin++)
1098 for (int col=0; col<nr; col++)
1099 ope.set(lin+nr,col) = -0.5*ms(lin,col) ;
1100 for (int lin=0; lin<nr; lin++)
1101 for (int col=0; col<nr; col++)
1102 ope.set(lin+nr,col+nr) = -md(lin,col) + 3*ms(lin,col) ;
1103 for (int lin=0; lin<nr; lin++)
1104 for (int col=0; col<nr; col++)
1105 ope.set(lin+nr,col+2*nr) = (2. - l_q*(l_q+1))*ms(lin,col) ;
1106 for (int lin=0; lin<nr; lin++)
1107 for (int col=0; col<nr; col++)
1108 ope.set(lin+2*nr,col) = 0.5*md(lin,col)/double(l_q+1)
1109 - 0.5*double(l_q+4)/double(l_q+1)*ms(lin,col) ;
1110 for (int lin=0; lin<nr; lin++)
1111 for (int col=0; col<nr; col++)
1112 ope.set(lin+2*nr,col+nr) = -2*ms(lin,col) ;
1113 for (int lin=0; lin<nr; lin++)
1114 for (int col=0; col<nr; col++)
1115 ope.set(lin+2*nr,col+2*nr) = -(l_q+2)*md(lin,col)
1116 + l_q*(l_q+2)*ms(lin,col) ;
1117 ope *= 1./alpha ;
1118 for (int col=0; col<3*nr; col++) {
1119 ope.set(ind0+nr-2, col) = 0 ;
1120 ope.set(ind0+nr-1, col) = 0 ;
1121 ope.set(ind1+nr-2, col) = 0 ;
1122 ope.set(ind1+nr-1, col) = 0 ;
1123 ope.set(ind2+nr-1, col) = 0 ;
1124 }
1125 for (int col=0; col<nr; col++) {
1126 ope.set(ind0+nr-1, col+ind0) = 1 ;
1127 ope.set(ind1+nr-1, col+ind1) = 1 ;
1128 ope.set(ind2+nr-1, col+ind2) = 1 ;
1129 }
1130 ope.set(ind0+nr-2, ind0+1) = 1 ;
1131 ope.set(ind1+nr-2, ind1+2) = 1 ;
1132
1133 ope.set_lu() ;
1134 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
1135 Matrice* pope = new Matrice(ope) ;
1136 par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
1137 mat_done.set(l_q) = 1 ;
1138 }
1139 } //End of case when a calculation is needed
1140 const Matrice& oper = (par_mat == 0x0 ? ope :
1141 par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
1142
1143 Tbl sec(3*nr) ;
1144 sec.set_etat_qcq() ;
1145 if (hnull) {
1146 for (int lin=0; lin<2*nr; lin++)
1147 sec.set(lin) = 0 ;
1148 for (int lin=0; lin<nr; lin++)
1149 sec.set(2*nr+lin) = (*source.get_spectral_va().c_cf)
1150 (lz, k, j, lin) ;
1151 }
1152 else {
1153 for (int lin=0; lin<nr; lin++)
1154 sec.set(lin) = (*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) ;
1155 for (int lin=0; lin<nr; lin++)
1156 sec.set(lin+nr) = -0.5*(*hoverr.get_spectral_va().c_cf)
1157 (lz, k, j, lin) ;
1158 if (bnull) {
1159 for (int lin=0; lin<nr; lin++)
1160 sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
1161 (*dhdr.get_spectral_va().c_cf)(lz, k, j, lin)
1162 + (l_q+2)*(*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) ) ;
1163 }
1164 else {
1165 for (int lin=0; lin<nr; lin++)
1166 sec.set(2*nr+lin) = -0.5/double(l_q+1)*(
1167 (*dhdr.get_spectral_va().c_cf)(lz, k, j, lin)
1168 + (l_q+2)*(*hoverr.get_spectral_va().c_cf)(lz, k, j, lin) )
1169 + (*source.get_spectral_va().c_cf)(lz, k, j, lin) ;
1170 }
1171 }
1172 sec.set(ind0+nr-2) = 0 ;
1173 sec.set(ind0+nr-1) = 0 ;
1174 sec.set(ind1+nr-1) = 0 ;
1175 sec.set(ind1+nr-2) = 0 ;
1176 sec.set(ind2+nr-1) = 0 ;
1177 Tbl sol = oper.inverse(sec) ;
1178 for (int i=0; i<nr; i++) {
1179 sol_part_hrr.set(lz, k, j, i) = sol(i) ;
1180 sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
1181 sol_part_w.set(lz, k, j, i) = sol(i+2*nr) ;
1182 }
1183 sec.annule_hard() ;
1184 sec.set(ind0+nr-2) = 1 ;
1185 sol = oper.inverse(sec) ;
1186 for (int i=0; i<nr; i++) {
1187 sol_hom1_hrr.set(lz, k, j, i) = sol(i) ;
1188 sol_hom1_eta.set(lz, k, j, i) = sol(i+nr) ;
1189 sol_hom1_w.set(lz, k, j, i) = sol(i+2*nr) ;
1190 }
1191 sec.set(ind0+nr-2) = 0 ;
1192 sec.set(ind1+nr-2) = 1 ;
1193 sol = oper.inverse(sec) ;
1194 for (int i=0; i<nr; i++) {
1195 sol_hom2_hrr.set(lz, k, j, i) = sol(i) ;
1196 sol_hom2_eta.set(lz, k, j, i) = sol(i+nr) ;
1197 sol_hom2_w.set(lz, k, j, i) = sol(i+2*nr) ;
1198 }
1199 }
1200 }
1201 }
1202 }
1203
1204 // Matching system...
1205
1206
1207
1208 int taille = 3*nz_bc ; // Un de moins que dans la version complete R3
1209 if (cedbc) taille-- ;
1210 Mtbl_cf& mhrr = *hrr.set_spectral_va().c_cf ;
1211 Mtbl_cf& meta = *tilde_eta.set_spectral_va().c_cf ;
1212 Mtbl_cf& mw = *ww.set_spectral_va().c_cf ;
1213 Tbl sec_membre(taille) ;
1214 Matrice systeme(taille, taille) ;
1215 int ligne ; int colonne ;
1216 Tbl pipo(1) ;
1217 const Tbl& hrrb = (cedbc ? pipo : par_bc->get_tbl_mod(1) );
1218 double chrr = (cedbc ? 0 : par_bc->get_tbl_mod()(4) ) ;
1219 double dhrr = (cedbc ? 0 : par_bc->get_tbl_mod()(5) ) ;
1220 double ceta = (cedbc ? 0 : par_bc->get_tbl_mod()(6) ) ;
1221 double deta = (cedbc ? 0 : par_bc->get_tbl_mod()(7) ) ;
1222 double cw = (cedbc ? 0 : par_bc->get_tbl_mod()(8) ) ;
1223 double dw = (cedbc ? 0 : par_bc->get_tbl_mod()(9) ) ;
1236
1237
1238 dhom1_hrr.dsdx() ; dhom1_eta.dsdx() ; dhom1_w.dsdx() ;
1239 dhom2_hrr.dsdx() ; dhom2_eta.dsdx() ; dhom2_w.dsdx() ;
1240 dhom3_hrr.dsdx() ; dhom3_eta.dsdx() ; dhom3_w.dsdx() ;
1241 dpart_hrr.dsdx() ; dpart_eta.dsdx() ; dpart_w.dsdx() ;
1242
1247 d2hom1_hrr.dsdx(); d2hom2_hrr.dsdx(); d2hom3_hrr.dsdx(); d2part_hrr.dsdx();
1248
1249
1250
1252 // Loop on l and m//
1254
1255 for (int k=0 ; k<np+1 ; k++)
1256 for (int j=0 ; j<nt ; j++) {
1257 base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
1258 if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
1259 ligne = 0 ;
1260 colonne = 0 ;
1261 systeme.annule_hard() ;
1262 sec_membre.annule_hard() ;
1263
1264 //First shell
1265 //Condition at x=-1;
1266 int nr = mgrid.get_nr(1);
1267 double alpha2 = mp_aff->get_alpha()[1] ;
1268 // A compatibility condition is set in order to be coherent with second order potentials definition
1269
1270 const Coord& rr = (*mp_aff).r ;
1271 Scalar rrr(*mp_aff); rrr = rr;
1272
1273 systeme.set(ligne, colonne)= (dhom1_w.val_in_bound_jk(1,j,k))/alpha2 + (0.5/rhor)*(double(l_q*(l_q+1)))*sol_hom1_w.val_in_bound_jk(1,j,k) - (1./rhor)*sol_hom1_eta.val_in_bound_jk(1,j,k) - (0.25/rhor)*sol_hom1_hrr.val_in_bound_jk(1,j,k);
1274
1275 systeme.set(ligne, colonne+1)= (dhom2_w.val_in_bound_jk(1,j,k))/alpha2 +(0.5/rhor)*(double(l_q*(l_q+1)))*sol_hom2_w.val_in_bound_jk(1,j,k) - (1./rhor)*sol_hom2_eta.val_in_bound_jk(1,j,k) - (0.25/rhor)*sol_hom2_hrr.val_in_bound_jk(1,j,k);
1276
1277 systeme.set(ligne, colonne+2)= (dhom3_w.val_in_bound_jk(1,j,k))/alpha2 +(0.5/rhor)*(double(l_q*(l_q+1)))*sol_hom3_w.val_in_bound_jk(1,j,k) - (1./rhor)*sol_hom3_eta.val_in_bound_jk(1,j,k) - (0.25/rhor)*sol_hom3_hrr.val_in_bound_jk(1,j,k) ;
1278
1279 sec_membre.set(ligne)= - (dpart_w.val_in_bound_jk(1,j,k)/alpha2 +(0.5/rhor)*(double(l_q*(l_q+1)))*sol_part_w.val_in_bound_jk(1,j,k) - (1./rhor)*sol_part_eta.val_in_bound_jk(1,j,k) - (0.25/rhor)*sol_part_hrr.val_in_bound_jk(1,j,k));
1280
1281 Mtbl_cf *Bbb = bb2.get_spectral_va().c_cf;
1282
1283 sec_membre.set(ligne) += (*Bbb).val_in_bound_jk(1,j,k);
1284
1285
1286 ligne++ ;
1287
1288 // A Robyn-type boundary condition is imposed on eta
1289
1290 systeme.set(ligne, colonne) =
1291 dir*sol_hom1_hrr.val_in_bound_jk(1, j, k) + neum*dhom1_hrr.val_in_bound_jk(1,j,k)/alpha2;
1292 systeme.set(ligne, colonne+1) =
1293 dir*sol_hom2_hrr.val_in_bound_jk(1, j, k) +neum*dhom2_hrr.val_in_bound_jk(1,j,k)/alpha2;
1294 systeme.set(ligne, colonne+2) =
1295 + dir*sol_hom3_hrr.val_in_bound_jk(1, j, k) + neum*dhom3_hrr.val_in_bound_jk(1,j,k)/alpha2;
1296
1297
1298 double blob = (*(bound_eta.get_spectral_va().c_cf)).val_in_bound_jk(1, j, k);
1299 sec_membre.set(ligne) += - sol_part_hrr.val_in_bound_jk(1, j, k) - neum*dpart_hrr.val_in_bound_jk(1,j,k)/alpha2 + blob ;
1300
1301 ligne++ ;
1302
1303
1304 //Conditions at x=1;
1305
1306 systeme.set(ligne, colonne) =
1307 sol_hom1_hrr.val_out_bound_jk(1, j, k) ;
1308 systeme.set(ligne, colonne+1) =
1309 sol_hom2_hrr.val_out_bound_jk(1, j, k) ;
1310 systeme.set(ligne, colonne+2) =
1311 sol_hom3_hrr.val_out_bound_jk(1, j, k) ;
1312
1313 sec_membre.set(ligne) -= sol_part_hrr.val_out_bound_jk(1, j, k) ;
1314 ligne++ ;
1315
1316 systeme.set(ligne, colonne) =
1317 sol_hom1_eta.val_out_bound_jk(1, j, k) ;
1318 systeme.set(ligne, colonne+1) =
1319 sol_hom2_eta.val_out_bound_jk(1, j, k) ;
1320 systeme.set(ligne, colonne+2) =
1321 sol_hom3_eta.val_out_bound_jk(1, j, k) ;
1322
1323 sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(1, j, k) ;
1324 ligne++ ;
1325
1326 systeme.set(ligne, colonne) =
1327 sol_hom1_w.val_out_bound_jk(1, j, k) ;
1328 systeme.set(ligne, colonne+1) =
1329 sol_hom2_w.val_out_bound_jk(1, j, k) ;
1330 systeme.set(ligne, colonne+2) =
1331 sol_hom3_w.val_out_bound_jk(1, j, k) ;
1332
1333 sec_membre.set(ligne) -= sol_part_w.val_out_bound_jk(1, j, k) ;
1334
1335 colonne += 3 ;
1336
1337 //shells
1338 for (int zone=2 ; zone<nz_bc ; zone++) {
1339 nr = mgrid.get_nr(zone) ;
1340 ligne -= 2 ;
1341
1342 //Condition at x = -1
1343 systeme.set(ligne, colonne) =
1344 - sol_hom1_hrr.val_in_bound_jk(zone, j, k) ;
1345 systeme.set(ligne, colonne+1) =
1346 - sol_hom2_hrr.val_in_bound_jk(zone, j, k) ;
1347 systeme.set(ligne, colonne+2) =
1348 - sol_hom3_hrr.val_in_bound_jk(zone, j, k) ;
1349
1350 sec_membre.set(ligne) += sol_part_hrr.val_in_bound_jk(zone, j, k) ;
1351 ligne++ ;
1352
1353 systeme.set(ligne, colonne) =
1354 - sol_hom1_eta.val_in_bound_jk(zone, j, k) ;
1355 systeme.set(ligne, colonne+1) =
1356 - sol_hom2_eta.val_in_bound_jk(zone, j, k) ;
1357 systeme.set(ligne, colonne+2) =
1358 - sol_hom3_eta.val_in_bound_jk(zone, j, k) ;
1359
1360 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(zone, j, k) ;
1361 ligne++ ;
1362
1363 systeme.set(ligne, colonne) =
1364 - sol_hom1_w.val_in_bound_jk(zone, j, k) ;
1365 systeme.set(ligne, colonne+1) =
1366 - sol_hom2_w.val_in_bound_jk(zone, j, k) ;
1367 systeme.set(ligne, colonne+2) =
1368 - sol_hom3_w.val_in_bound_jk(zone, j, k) ;
1369
1370 sec_membre.set(ligne) += sol_part_w.val_in_bound_jk(zone, j, k) ;
1371 ligne++ ;
1372
1373 // Condition at x=1
1374 systeme.set(ligne, colonne) =
1375 sol_hom1_hrr.val_out_bound_jk(zone, j, k) ;
1376 systeme.set(ligne, colonne+1) =
1377 sol_hom2_hrr.val_out_bound_jk(zone, j, k) ;
1378 systeme.set(ligne, colonne+2) =
1379 sol_hom3_hrr.val_out_bound_jk(zone, j, k) ;
1380
1381 sec_membre.set(ligne) -= sol_part_hrr.val_out_bound_jk(zone, j, k) ;
1382 ligne++ ;
1383
1384 systeme.set(ligne, colonne) =
1385 sol_hom1_eta.val_out_bound_jk(zone, j, k) ;
1386 systeme.set(ligne, colonne+1) =
1387 sol_hom2_eta.val_out_bound_jk(zone, j, k) ;
1388 systeme.set(ligne, colonne+2) =
1389 sol_hom3_eta.val_out_bound_jk(zone, j, k) ;
1390
1391 sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone, j, k) ;
1392 ligne++ ;
1393
1394 systeme.set(ligne, colonne) =
1395 sol_hom1_w.val_out_bound_jk(zone, j, k) ;
1396 systeme.set(ligne, colonne+1) =
1397 sol_hom2_w.val_out_bound_jk(zone, j, k) ;
1398 systeme.set(ligne, colonne+2) =
1399 sol_hom3_w.val_out_bound_jk(zone, j, k) ;
1400
1401 sec_membre.set(ligne) -= sol_part_w.val_out_bound_jk(zone, j, k) ;
1402
1403 colonne += 3 ;
1404 }
1405
1406 //Last domain
1407 nr = mgrid.get_nr(nz_bc) ;
1408 double alpha = mp_aff->get_alpha()[nz_bc] ;
1409 ligne -= 2 ;
1410
1411 //Condition at x = -1
1412 systeme.set(ligne, colonne) =
1413 - sol_hom1_hrr.val_in_bound_jk(nz_bc, j, k) ;
1414 systeme.set(ligne, colonne+1) =
1415 - sol_hom2_hrr.val_in_bound_jk(nz_bc, j, k) ;
1416 if (!cedbc) systeme.set(ligne, colonne+2) =
1417 - sol_hom3_hrr.val_in_bound_jk(nz_bc, j, k) ;
1418
1419 sec_membre.set(ligne) += sol_part_hrr.val_in_bound_jk(nz_bc, j, k) ;
1420 ligne++ ;
1421
1422 systeme.set(ligne, colonne) =
1423 - sol_hom1_eta.val_in_bound_jk(nz_bc, j, k) ;
1424 systeme.set(ligne, colonne+1) =
1425 - sol_hom2_eta.val_in_bound_jk(nz_bc, j, k) ;
1426 if (!cedbc) systeme.set(ligne, colonne+2) =
1427 - sol_hom3_eta.val_in_bound_jk(nz_bc, j, k) ;
1428
1429 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(nz_bc, j, k) ;
1430 ligne++ ;
1431
1432 systeme.set(ligne, colonne) =
1433 - sol_hom1_w.val_in_bound_jk(nz_bc, j, k) ;
1434 systeme.set(ligne, colonne+1) =
1435 - sol_hom2_w.val_in_bound_jk(nz_bc, j, k) ;
1436 if (!cedbc) systeme.set(ligne, colonne+2) =
1437 - sol_hom3_w.val_in_bound_jk(nz_bc, j, k) ;
1438
1439 sec_membre.set(ligne) += sol_part_w.val_in_bound_jk(nz_bc, j, k) ;
1440 ligne++ ;
1441
1442 if (!cedbc) {//Special condition at x=1
1443 systeme.set(ligne, colonne) =
1444 chrr*sol_hom1_hrr.val_out_bound_jk(nz_bc, j, k)
1445 + dhrr*dhom1_hrr.val_out_bound_jk(nz_bc, j, k) / alpha
1446 + ceta*sol_hom1_eta.val_out_bound_jk(nz_bc, j, k)
1447 + deta*dhom1_eta.val_out_bound_jk(nz_bc, j, k) / alpha
1448 + cw*sol_hom1_w.val_out_bound_jk(nz_bc, j, k)
1449 + dw*dhom1_w.val_out_bound_jk(nz_bc, j, k) / alpha ;
1450 systeme.set(ligne, colonne+1) =
1451 chrr*sol_hom2_hrr.val_out_bound_jk(nz_bc, j, k)
1452 + dhrr*dhom2_hrr.val_out_bound_jk(nz_bc, j, k) / alpha
1453 + ceta*sol_hom2_eta.val_out_bound_jk(nz_bc, j, k)
1454 + deta*dhom2_eta.val_out_bound_jk(nz_bc, j, k) / alpha
1455 + cw*sol_hom2_w.val_out_bound_jk(nz_bc, j, k)
1456 + dw*dhom2_w.val_out_bound_jk(nz_bc, j, k) / alpha ;
1457 systeme.set(ligne, colonne+2) =
1458 chrr*sol_hom3_hrr.val_out_bound_jk(nz_bc, j, k)
1459 + dhrr*dhom3_hrr.val_out_bound_jk(nz_bc, j, k) / alpha
1460 + ceta*sol_hom3_eta.val_out_bound_jk(nz_bc, j, k)
1461 + deta*dhom3_eta.val_out_bound_jk(nz_bc, j, k) / alpha
1462 + cw*sol_hom3_w.val_out_bound_jk(nz_bc, j, k)
1463 + dw*dhom3_w.val_out_bound_jk(nz_bc, j, k) / alpha ;
1464
1465 sec_membre.set(ligne) -= chrr*sol_part_hrr.val_out_bound_jk(nz_bc, j, k)
1466 + dhrr*dpart_hrr.val_out_bound_jk(nz_bc, j, k)/alpha
1467 + ceta*sol_part_eta.val_out_bound_jk(nz_bc, j, k)
1468 + deta*dpart_eta.val_out_bound_jk(nz_bc, j, k)/alpha
1469 + cw*sol_part_w.val_out_bound_jk(nz_bc, j, k)
1470 + dw*dpart_w.val_out_bound_jk(nz_bc, j, k)/alpha
1471 - hrrb(k, j) ;
1472 }
1473
1474
1476
1477
1478
1479 //System inversion: Solution of the system giving the coefficients for the homogeneous
1480 // solutions
1481 //-------------------------------------------------------------------
1482
1483 systeme.set_lu() ;
1484 Tbl facteur = systeme.inverse(sec_membre) ;
1485 int conte = 0 ;
1486
1487
1488 // everything is put to the right place,
1489 //---------------------------------------
1490 nr = mgrid.get_nr(0) ; //nucleus
1491 for (int i=0 ; i<nr ; i++) {
1492 mhrr.set(0, k, j, i) = 0 ;
1493 meta.set(0, k, j, i) = 0 ;
1494 mw.set(0, k, j, i) = 0 ;
1495 }
1496 for (int zone=1 ; zone<=n_shell ; zone++) { //shells
1497 nr = mgrid.get_nr(zone) ;
1498 for (int i=0 ; i<nr ; i++) {
1499 mhrr.set(zone, k, j, i) = sol_part_hrr(zone, k, j, i)
1501 + facteur(conte+1)*sol_hom2_hrr(zone, k, j, i)
1502 + facteur(conte+2)*sol_hom3_hrr(zone, k, j, i) ;
1503
1504 meta.set(zone, k, j, i) = sol_part_eta(zone, k, j, i)
1506 + facteur(conte+1)*sol_hom2_eta(zone, k, j, i)
1507 + facteur(conte+2)*sol_hom3_eta(zone, k, j, i) ;
1508
1509 mw.set(zone, k, j, i) = sol_part_w(zone, k, j, i)
1510 + facteur(conte)*sol_hom1_w(zone, k, j, i)
1511 + facteur(conte+1)*sol_hom2_w(zone, k, j, i)
1512 + facteur(conte+2)*sol_hom3_w(zone, k, j, i) ;
1513 }
1514 conte+=3 ;
1515 }
1516 if (cedbc) {
1517 nr = mgrid.get_nr(nzm1) ; //compactified external domain
1518 for (int i=0 ; i<nr ; i++) {
1519 mhrr.set(nzm1, k, j, i) = sol_part_hrr(nzm1, k, j, i)
1521 + facteur(conte+1)*sol_hom2_hrr(nzm1, k, j, i) ;
1522
1523 meta.set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
1525 + facteur(conte+1)*sol_hom2_eta(nzm1, k, j, i) ;
1526
1527 mw.set(nzm1, k, j, i) = sol_part_w(nzm1, k, j, i)
1528 + facteur(conte)*sol_hom1_w(nzm1, k, j, i)
1529 + facteur(conte+1)*sol_hom2_w(nzm1, k, j, i) ;
1530 }
1531 }
1532 } // End of nullite_plm
1533 } //End of loop on theta
1534
1535
1536 if (hrr.set_spectral_va().c != 0x0)
1537 delete hrr.set_spectral_va().c;
1538 hrr.set_spectral_va().c = 0x0 ;
1539 hrr.set_spectral_va().ylm_i() ;
1540
1541 if (tilde_eta.set_spectral_va().c != 0x0)
1542 delete tilde_eta.set_spectral_va().c ;
1543 tilde_eta.set_spectral_va().c = 0x0 ;
1544 tilde_eta.set_spectral_va().ylm_i() ;
1545
1546 if (ww.set_spectral_va().c != 0x0)
1547 delete ww.set_spectral_va().c ;
1548 ww.set_spectral_va().c = 0x0 ;
1549 ww.set_spectral_va().ylm_i() ;
1550
1551
1552
1553// // On doit resoudre séparément les cas l=0 et l=1; notamment dansces cas le systeme electrique se resout a deux equations
1554// // au lieu de 3.
1555
1556
1557}
1558
1559void Sym_tensor_trans::sol_Dirac_l01_2(const Scalar& hh, Scalar& hrr, Scalar& tilde_eta,
1560 Param* par_mat) {
1561
1562
1563
1564 // void Sym_tensor_trans::sol_Dirac_tilde_B2(const Scalar& tilde_b, const Scalar& hh,
1565 // Scalar& hrr, Scalar& tilde_eta, Scalar& ww, Scalar bound_eta,double dir, double neum, double rhor, Param* par_bc, Param* par_mat) {
1566
1567
1568 const Map_af* mp_aff = dynamic_cast<const Map_af*>(mp) ;
1569 assert(mp_aff != 0x0) ; //Only affine mapping for the moment
1570
1571 const Mg3d& mgrid = *mp_aff->get_mg() ;
1572 int nz = mgrid.get_nzone() ;
1573 assert(mgrid.get_type_r(0) == RARE) ;
1574 assert(mgrid.get_type_r(nz-1) == UNSURR) ;
1575
1576 if (hh.get_etat() == ETATZERO) {
1577 hrr.annule_hard() ;
1578 tilde_eta.annule_hard() ;
1579 return ;
1580 }
1581
1582 int nt = mgrid.get_nt(0) ;
1583 int np = mgrid.get_np(0) ;
1584
1585 Scalar source = hh ;
1586 source.div_r_dzpuis(2) ;
1587 Scalar source_coq = hh ;
1588 source.set_spectral_va().ylm() ;
1589 source_coq.set_spectral_va().ylm() ;
1590 Base_val base = source.get_spectral_base() ;
1591 base.mult_x() ;
1592 int lmax = base.give_lmax(mgrid, 0) + 1;
1593
1594 assert (hrr.get_spectral_base() == base) ;
1595 assert (tilde_eta.get_spectral_base() == base) ;
1596 assert (hrr.get_spectral_va().c_cf != 0x0) ;
1597 assert (tilde_eta.get_spectral_va().c_cf != 0x0) ;
1598
1599 Mtbl_cf sol_part_hrr(mgrid, base) ; sol_part_hrr.annule_hard() ;
1600 Mtbl_cf sol_part_eta(mgrid, base) ; sol_part_eta.annule_hard() ;
1601 Mtbl_cf sol_hom1_hrr(mgrid, base) ; sol_hom1_hrr.annule_hard() ;
1602 Mtbl_cf sol_hom1_eta(mgrid, base) ; sol_hom1_eta.annule_hard() ;
1603 Mtbl_cf sol_hom2_hrr(mgrid, base) ; sol_hom2_hrr.annule_hard() ;
1604 Mtbl_cf sol_hom2_eta(mgrid, base) ; sol_hom2_eta.annule_hard() ;
1605
1606 bool need_calculation = true ;
1607 if (par_mat != 0x0)
1608 if (par_mat->get_n_matrice_mod() > 0)
1609 if (&par_mat->get_matrice_mod(0) != 0x0) need_calculation = false ;
1610
1611 int l_q, m_q, base_r ;
1612 Itbl mat_done(lmax) ;
1613
1614 //---------------
1615 //-- NUCLEUS ---
1616 //---------------
1617 {int lz = 0 ;
1618 int nr = mgrid.get_nr(lz) ;
1619 double alpha = mp_aff->get_alpha()[lz] ;
1620 Matrice ope(2*nr, 2*nr) ;
1621 if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
1622
1623 for (int k=0 ; k<np+1 ; k++) {
1624 for (int j=0 ; j<nt ; j++) {
1625 // quantic numbers and spectral bases
1626 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
1627 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q < 2)) {
1628 if (need_calculation) {
1629 ope.set_etat_qcq() ;
1630 Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
1631 Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
1632
1633 for (int lin=0; lin<nr; lin++)
1634 for (int col=0; col<nr; col++)
1635 ope.set(lin,col) = md(lin,col) + 3*ms(lin,col) ;
1636 for (int lin=0; lin<nr; lin++)
1637 for (int col=0; col<nr; col++)
1638 ope.set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
1639 for (int lin=0; lin<nr; lin++)
1640 for (int col=0; col<nr; col++)
1641 ope.set(lin+nr,col) = -0.5*ms(lin,col) ;
1642 for (int lin=0; lin<nr; lin++)
1643 for (int col=0; col<nr; col++)
1644 ope.set(lin+nr,col+nr) = md(lin,col) + 3*ms(lin, col);
1645
1646 ope *= 1./alpha ;
1647 for (int col=0; col<2*nr; col++) {
1648 ope.set(nr-1, col) = 0 ;
1649 ope.set(2*nr-1, col) = 0 ;
1650 }
1651 int pari = 1 ;
1652 if (base_r == R_CHEBP) {
1653 for (int col=0; col<nr; col++) {
1654 ope.set(nr-1, col) = pari ;
1655 ope.set(2*nr-1, col+nr) = pari ;
1656 pari = - pari ;
1657 }
1658 }
1659 else { //In the odd case, the last coefficient must be zero!
1660 ope.set(nr-1, nr-1) = 1 ;
1661 ope.set(2*nr-1, 2*nr-1) = 1 ;
1662 }
1663
1664 ope.set_lu() ;
1665 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
1666 Matrice* pope = new Matrice(ope) ;
1667 par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
1668 mat_done.set(l_q) = 1 ;
1669 }
1670 } //End of case when a calculation is needed
1671
1672 const Matrice& oper = (par_mat == 0x0 ? ope :
1673 par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
1674 Tbl sec(2*nr) ;
1675 sec.set_etat_qcq() ;
1676 for (int lin=0; lin<nr; lin++)
1677 sec.set(lin) = (*source.get_spectral_va().c_cf)(lz, k, j, lin) ;
1678 for (int lin=0; lin<nr; lin++)
1679 sec.set(nr+lin) = -0.5*(*source.get_spectral_va().c_cf)
1680 (lz, k, j, lin) ;
1681 sec.set(nr-1) = 0 ;
1682 if (base_r == R_CHEBP) {
1683 double h0 = 0 ; //In the l=0 case: 3*hrr(r=0) = h(r=0)
1684 int pari = 1 ;
1685 for (int col=0; col<nr; col++) {
1686 h0 += pari*
1687 (*source_coq.get_spectral_va().c_cf)(lz, k, j, col) ;
1688 pari = - pari ;
1689 }
1690 sec.set(nr-1) = h0 / 3. ;
1691 }
1692 sec.set(2*nr-1) = 0 ;
1693 Tbl sol = oper.inverse(sec) ;
1694 for (int i=0; i<nr; i++) {
1695 sol_part_hrr.set(lz, k, j, i) = sol(i) ;
1696 sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
1697 }
1698 sec.annule_hard() ;
1699 }
1700 }
1701 }
1702 }
1703
1704 //-------------
1705 // -- Shells --
1706 //-------------
1707
1708 for (int lz=1; lz<nz-1; lz++) {
1709 if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
1710 int nr = mgrid.get_nr(lz) ;
1711 int ind0 = 0 ;
1712 int ind1 = nr ;
1713 assert(mgrid.get_nt(lz) == nt) ;
1714 assert(mgrid.get_np(lz) == np) ;
1715 double alpha = mp_aff->get_alpha()[lz] ;
1716 double ech = mp_aff->get_beta()[lz] / alpha ;
1717 Matrice ope(2*nr, 2*nr) ;
1718
1719 for (int k=0 ; k<np+1 ; k++) {
1720 for (int j=0 ; j<nt ; j++) {
1721 // quantic numbers and spectral bases
1722 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
1723 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q < 2)) {
1724 if (need_calculation) {
1725 ope.set_etat_qcq() ;
1726 Diff_xdsdx oxd(base_r, nr) ; const Matrice& mxd = oxd.get_matrice() ;
1727 Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
1728 Diff_id oid(base_r, nr) ; const Matrice& mid = oid.get_matrice() ;
1729
1730 for (int lin=0; lin<nr; lin++)
1731 for (int col=0; col<nr; col++)
1732 ope.set(lin,col) = mxd(lin,col) + ech*md(lin,col)
1733 + 3*mid(lin,col) ;
1734 for (int lin=0; lin<nr; lin++)
1735 for (int col=0; col<nr; col++)
1736 ope.set(lin,col+nr) = -l_q*(l_q+1)*mid(lin,col) ;
1737 for (int lin=0; lin<nr; lin++)
1738 for (int col=0; col<nr; col++)
1739 ope.set(lin+nr,col) = -0.5*mid(lin,col) ;
1740 for (int lin=0; lin<nr; lin++)
1741 for (int col=0; col<nr; col++)
1742 ope.set(lin+nr,col+nr) = mxd(lin,col) + ech*md(lin,col)
1743 + 3*mid(lin, col) ;
1744
1745 for (int col=0; col<2*nr; col++) {
1746 ope.set(ind0+nr-1, col) = 0 ;
1747 ope.set(ind1+nr-1, col) = 0 ;
1748 }
1749 ope.set(ind0+nr-1, ind0) = 1 ;
1750 ope.set(ind1+nr-1, ind1) = 1 ;
1751
1752 ope.set_lu() ;
1753 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
1754 Matrice* pope = new Matrice(ope) ;
1755 par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
1756 mat_done.set(l_q) = 1 ;
1757 }
1758 } //End of case when a calculation is needed
1759 const Matrice& oper = (par_mat == 0x0 ? ope :
1760 par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
1761 Tbl sec(2*nr) ;
1762 sec.set_etat_qcq() ;
1763 for (int lin=0; lin<nr; lin++)
1764 sec.set(lin) = (*source_coq.get_spectral_va().c_cf)
1765 (lz, k, j, lin) ;
1766 for (int lin=0; lin<nr; lin++)
1767 sec.set(nr+lin) = -0.5*(*source_coq.get_spectral_va().c_cf)
1768 (lz, k, j, lin) ;
1769 sec.set(ind0+nr-1) = 0 ;
1770 sec.set(ind1+nr-1) = 0 ;
1771 Tbl sol = oper.inverse(sec) ;
1772
1773 for (int i=0; i<nr; i++) {
1774 sol_part_hrr.set(lz, k, j, i) = sol(i) ;
1775 sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
1776 }
1777 sec.annule_hard() ;
1778 sec.set(ind0+nr-1) = 1 ;
1779 sol = oper.inverse(sec) ;
1780 for (int i=0; i<nr; i++) {
1781 sol_hom1_hrr.set(lz, k, j, i) = sol(i) ;
1782 sol_hom1_eta.set(lz, k, j, i) = sol(i+nr) ;
1783 }
1784 sec.set(ind0+nr-1) = 0 ;
1785 sec.set(ind1+nr-1) = 1 ;
1786 sol = oper.inverse(sec) ;
1787 for (int i=0; i<nr; i++) {
1788 sol_hom2_hrr.set(lz, k, j, i) = sol(i) ;
1789 sol_hom2_eta.set(lz, k, j, i) = sol(i+nr) ;
1790 }
1791 }
1792 }
1793 }
1794 }
1795
1796 //------------------------------
1797 // Compactified external domain
1798 //------------------------------
1799 {int lz = nz-1 ;
1800 if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
1801 int nr = mgrid.get_nr(lz) ;
1802 int ind0 = 0 ;
1803 int ind1 = nr ;
1804 assert(mgrid.get_nt(lz) == nt) ;
1805 assert(mgrid.get_np(lz) == np) ;
1806 double alpha = mp_aff->get_alpha()[lz] ;
1807 Matrice ope(2*nr, 2*nr) ;
1808
1809 for (int k=0 ; k<np+1 ; k++) {
1810 for (int j=0 ; j<nt ; j++) {
1811 // quantic numbers and spectral bases
1812 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
1813 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q < 2)) {
1814 if (need_calculation) {
1815 ope.set_etat_qcq() ;
1816 Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
1817 Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
1818
1819 for (int lin=0; lin<nr; lin++)
1820 for (int col=0; col<nr; col++)
1821 ope.set(lin,col) = - md(lin,col) + 3*ms(lin,col) ;
1822 for (int lin=0; lin<nr; lin++)
1823 for (int col=0; col<nr; col++)
1824 ope.set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
1825 for (int lin=0; lin<nr; lin++)
1826 for (int col=0; col<nr; col++)
1827 ope.set(lin+nr,col) = -0.5*ms(lin,col) ;
1828 for (int lin=0; lin<nr; lin++)
1829 for (int col=0; col<nr; col++)
1830 ope.set(lin+nr,col+nr) = -md(lin,col) + 3*ms(lin, col) ;
1831
1832 ope *= 1./alpha ;
1833 for (int col=0; col<2*nr; col++) {
1834 ope.set(ind0+nr-2, col) = 0 ;
1835 ope.set(ind0+nr-1, col) = 0 ;
1836 ope.set(ind1+nr-2, col) = 0 ;
1837 ope.set(ind1+nr-1, col) = 0 ;
1838 }
1839 for (int col=0; col<nr; col++) {
1840 ope.set(ind0+nr-1, col+ind0) = 1 ;
1841 ope.set(ind1+nr-1, col+ind1) = 1 ;
1842 }
1843 ope.set(ind0+nr-2, ind0+1) = 1 ;
1844 ope.set(ind1+nr-2, ind1+1) = 1 ;
1845
1846 ope.set_lu() ;
1847 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
1848 Matrice* pope = new Matrice(ope) ;
1849 par_mat->add_matrice_mod(*pope, lz*lmax + l_q) ;
1850 mat_done.set(l_q) = 1 ;
1851 }
1852 } //End of case when a calculation is needed
1853 const Matrice& oper = (par_mat == 0x0 ? ope :
1854 par_mat->get_matrice_mod(lz*lmax + l_q) ) ;
1855 Tbl sec(2*nr) ;
1856 sec.set_etat_qcq() ;
1857 for (int lin=0; lin<nr; lin++)
1858 sec.set(lin) = (*source.get_spectral_va().c_cf)
1859 (lz, k, j, lin) ;
1860 for (int lin=0; lin<nr; lin++)
1861 sec.set(nr+lin) = -0.5*(*source.get_spectral_va().c_cf)
1862 (lz, k, j, lin) ;
1863 sec.set(ind0+nr-2) = 0 ;
1864 sec.set(ind0+nr-1) = 0 ;
1865 sec.set(ind1+nr-2) = 0 ;
1866 sec.set(ind1+nr-1) = 0 ;
1867 Tbl sol = oper.inverse(sec) ;
1868 for (int i=0; i<nr; i++) {
1869 sol_part_hrr.set(lz, k, j, i) = sol(i) ;
1870 sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
1871 }
1872 sec.annule_hard() ;
1873 sec.set(ind0+nr-2) = 1 ;
1874 sol = oper.inverse(sec) ;
1875 for (int i=0; i<nr; i++) {
1876 sol_hom1_hrr.set(lz, k, j, i) = sol(i) ;
1877 sol_hom1_eta.set(lz, k, j, i) = sol(i+nr) ;
1878 }
1879 sec.set(ind0+nr-2) = 0 ;
1880 sec.set(ind1+nr-2) = 1 ;
1881 sol = oper.inverse(sec) ;
1882 for (int i=0; i<nr; i++) {
1883 sol_hom2_hrr.set(lz, k, j, i) = sol(i) ;
1884 sol_hom2_eta.set(lz, k, j, i) = sol(i+nr) ;
1885 }
1886 }
1887 }
1888 }
1889 }
1890
1891 int taille = 2*(nz-1) ;
1892 Mtbl_cf& mhrr = *hrr.set_spectral_va().c_cf ;
1893 Mtbl_cf& meta = *tilde_eta.set_spectral_va().c_cf ;
1894
1895 Tbl sec_membre(taille) ;
1896 Matrice systeme(taille, taille) ;
1897 int ligne ; int colonne ;
1898
1899 // Loop on l and m
1900 //----------------
1901 for (int k=0 ; k<np+1 ; k++)
1902 for (int j=0 ; j<nt ; j++) {
1903 base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
1904 if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q < 2)) {
1905 ligne = 0 ;
1906 colonne = 0 ;
1907 systeme.annule_hard() ;
1908 sec_membre.annule_hard() ;
1909
1910 //Nucleus
1911 int nr = mgrid.get_nr(0) ;
1912
1913 sec_membre.set(ligne) = -sol_part_hrr.val_out_bound_jk(0, j, k) ;
1914 ligne++ ;
1915
1916 sec_membre.set(ligne) = -sol_part_eta.val_out_bound_jk(0, j, k) ;
1917
1918 //shells
1919 for (int zone=1 ; zone<nz-1 ; zone++) {
1920 nr = mgrid.get_nr(zone) ;
1921 ligne-- ;
1922
1923 //Condition at x = -1
1924 systeme.set(ligne, colonne) =
1925 - sol_hom1_hrr.val_in_bound_jk(zone, j, k) ;
1926 systeme.set(ligne, colonne+1) =
1927 - sol_hom2_hrr.val_in_bound_jk(zone, j, k) ;
1928
1929 sec_membre.set(ligne) += sol_part_hrr.val_in_bound_jk(zone, j, k) ;
1930 ligne++ ;
1931
1932 systeme.set(ligne, colonne) =
1933 - sol_hom1_eta.val_in_bound_jk(zone, j, k) ;
1934 systeme.set(ligne, colonne+1) =
1935 - sol_hom2_eta.val_in_bound_jk(zone, j, k) ;
1936
1937 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(zone, j, k) ;
1938 ligne++ ;
1939
1940 // Condition at x=1
1941 systeme.set(ligne, colonne) =
1942 sol_hom1_hrr.val_out_bound_jk(zone, j, k) ;
1943 systeme.set(ligne, colonne+1) =
1944 sol_hom2_hrr.val_out_bound_jk(zone, j, k) ;
1945
1946 sec_membre.set(ligne) -= sol_part_hrr.val_out_bound_jk(zone, j, k) ;
1947 ligne++ ;
1948
1949 systeme.set(ligne, colonne) =
1950 sol_hom1_eta.val_out_bound_jk(zone, j, k) ;
1951 systeme.set(ligne, colonne+1) =
1952 sol_hom2_eta.val_out_bound_jk(zone, j, k) ;
1953
1954 sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone, j, k) ;
1955
1956 colonne += 2 ;
1957 }
1958
1959 //Compactified external domain
1960 nr = mgrid.get_nr(nz-1) ;
1961
1962 ligne-- ;
1963
1964 systeme.set(ligne, colonne) =
1965 - sol_hom1_hrr.val_in_bound_jk(nz-1, j, k) ;
1966 systeme.set(ligne, colonne+1) =
1967 - sol_hom2_hrr.val_in_bound_jk(nz-1, j, k) ;
1968
1969 sec_membre.set(ligne) += sol_part_hrr.val_in_bound_jk(nz-1, j, k) ;
1970 ligne++ ;
1971
1972 systeme.set(ligne, colonne) =
1973 - sol_hom1_eta.val_in_bound_jk(nz-1, j, k) ;
1974 systeme.set(ligne, colonne+1) =
1975 - sol_hom2_eta.val_in_bound_jk(nz-1, j, k) ;
1976
1977 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(nz-1, j, k) ;
1978
1979 // Solution of the system giving the coefficients for the homogeneous
1980 // solutions
1981 //-------------------------------------------------------------------
1982 systeme.set_lu() ;
1983 Tbl facteur = systeme.inverse(sec_membre) ;
1984 int conte = 0 ;
1985
1986 // everything is put to the right place...
1987 //----------------------------------------
1988 nr = mgrid.get_nr(0) ; //nucleus
1989 for (int i=0 ; i<nr ; i++) {
1990 mhrr.set(0, k, j, i) = sol_part_hrr(0, k, j, i) ;
1991 meta.set(0, k, j, i) = sol_part_eta(0, k, j, i) ;
1992 }
1993 for (int zone=1 ; zone<nz-1 ; zone++) { //shells
1994 nr = mgrid.get_nr(zone) ;
1995 for (int i=0 ; i<nr ; i++) {
1996 mhrr.set(zone, k, j, i) = sol_part_hrr(zone, k, j, i)
1997 + facteur(conte)*sol_hom1_hrr(zone, k, j, i)
1998 + facteur(conte+1)*sol_hom2_hrr(zone, k, j, i) ;
1999
2000 meta.set(zone, k, j, i) = sol_part_eta(zone, k, j, i)
2001 + facteur(conte)*sol_hom1_eta(zone, k, j, i)
2002 + facteur(conte+1)*sol_hom2_eta(zone, k, j, i) ;
2003 }
2004 conte+=2 ;
2005 }
2006 nr = mgrid.get_nr(nz-1) ; //compactified external domain
2007 for (int i=0 ; i<nr ; i++) {
2008 mhrr.set(nz-1, k, j, i) = sol_part_hrr(nz-1, k, j, i)
2009 + facteur(conte)*sol_hom1_hrr(nz-1, k, j, i)
2010 + facteur(conte+1)*sol_hom2_hrr(nz-1, k, j, i) ;
2011
2012 meta.set(nz-1, k, j, i) = sol_part_eta(nz-1, k, j, i)
2013 + facteur(conte)*sol_hom1_eta(nz-1, k, j, i)
2014 + facteur(conte+1)*sol_hom2_eta(nz-1, k, j, i) ;
2015 }
2016
2017 } // End of nullite_plm
2018 } //End of loop on theta
2019}
2020
2021}
Bases of the spectral expansions.
Definition base_val.h:322
void mult_x()
The basis is transformed as with a multiplication by .
int give_lmax(const Mg3d &mgrid, int lz) const
Returns the highest multipole for a given grid.
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:129
Class for the elementary differential operator Identity (see the base class Diff ).
Definition diff.h:210
Class for the elementary differential operator division by (see the base class Diff ).
Definition diff.h:329
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:409
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Basic integer array class.
Definition itbl.h:122
Affine radial mapping.
Definition map.h:2027
const double * get_beta() const
Returns the pointer on the array beta.
Definition map_af.C:481
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
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition mtbl_cf.C:312
Parameter storage.
Definition param.h:125
int get_n_matrice_mod() const
Returns the number of modifiable Matrice 's addresses in the list.
Definition param.C:859
void add_matrice_mod(Matrice &ti, int position=0)
Adds the address of a new modifiable Matrice to the list.
Definition param.C:866
Matrice & get_matrice_mod(int position=0) const
Returns the reference of a modifiable Matrice stored in the list.
Definition param.C:911
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
void div_r_dzpuis(int ced_mult_r)
Division by r everywhere but with the output flag dzpuis set to ced_mult_r .
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition scalar.C:873
const Scalar & dsdr() const
Returns of *this .
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition scalar.h:601
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va
Definition scalar.h:1294
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
Definition scalar.C:797
void sol_Dirac_BC2(const Scalar &bb, const Scalar &cc, const Scalar &hh, Scalar &hrr, Scalar &tilde_eta, Scalar &ww, Scalar bound_eta, double dir, double neum, double rhor, Param *par_bc, Param *par_mat)
Same resolution as sol_Dirac_tilde_B, but with inner boundary conditions added.
void sol_Dirac_Abound(const Scalar &aaa, Scalar &tilde_mu, Scalar &x_new, Scalar bound_mu, const Param *par_bc)
Same resolution as sol_Dirac_A, but with inner boundary conditions added.
void sol_Dirac_l01(const Scalar &hh, Scalar &hrr, Scalar &tilde_eta, Param *par_mat) const
Solves the same system as Sym_tensor_trans::sol_Dirac_tilde_B but only for .
Basic array class.
Definition tbl.h:161
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
Definition valeur.C:712
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 ylm_i()
Inverse of ylm()
#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