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