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