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