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