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