LORENE
vector_divfree_A.C
1/*
2 * Methods to impose the Dirac gauge: divergence-free condition.
3 *
4 * (see file sym_tensor.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2006 Jerome Novak
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28
29
30/*
31 * $Id: vector_divfree_A.C,v 1.7 2016/12/05 16:18:18 j_novak Exp $
32 * $Log: vector_divfree_A.C,v $
33 * Revision 1.7 2016/12/05 16:18:18 j_novak
34 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35 *
36 * Revision 1.6 2014/10/13 08:53:44 j_novak
37 * Lorene classes and functions now belong to the namespace Lorene.
38 *
39 * Revision 1.5 2014/10/06 15:13:20 j_novak
40 * Modified #include directives to use c++ syntax.
41 *
42 * Revision 1.4 2013/06/05 15:10:43 j_novak
43 * Suppression of FINJAC sampling in r. This Jacobi(0,2) base is now
44 * available by setting colloc_r to BASE_JAC02 in the Mg3d constructor.
45 *
46 * Revision 1.3 2009/10/23 13:18:46 j_novak
47 * Minor modifications
48 *
49 * Revision 1.2 2008/08/27 10:55:15 jl_cornou
50 * Added the case of one zone, which is a nucleus for BC
51 *
52 * Revision 1.1 2008/08/27 09:01:27 jl_cornou
53 * Methods for solving Dirac systems for divergence free vectors
54 *
55 *
56 *
57 * $Header: /cvsroot/Lorene/C++/Source/Tensor/vector_divfree_A.C,v 1.7 2016/12/05 16:18:18 j_novak Exp $
58 *
59 */
60
61
62// C headers
63#include <cstdlib>
64#include <cassert>
65#include <cmath>
66
67// Lorene headers
68#include "metric.h"
69#include "diff.h"
70#include "proto.h"
71#include "param.h"
72
73//----------------------------------------------------------------------------------
74//
75// sol_Dirac_A
76//
77//----------------------------------------------------------------------------------
78
79namespace Lorene {
80void Vector_divfree::sol_Dirac_A(const Scalar& aaa, Scalar& tilde_vr, Scalar& tilde_eta,
81 const Param* par_bc) const {
82
83 const Map_af* mp_aff = dynamic_cast<const Map_af*>(mp) ;
84 assert(mp_aff != 0x0) ; //Only affine mapping for the moment
85
86 const Mg3d& mgrid = *mp_aff->get_mg() ;
87 assert(mgrid.get_type_r(0) == RARE) ;
88 if (aaa.get_etat() == ETATZERO) {
89 tilde_vr = 0 ;
90 tilde_eta = 0 ;
91 return ;
92 }
93 assert(aaa.get_etat() != ETATNONDEF) ;
94 int nz = mgrid.get_nzone() ;
95 int nzm1 = nz - 1 ;
96 bool ced = (mgrid.get_type_r(nzm1) == UNSURR) ;
97 int n_shell = ced ? nz-2 : nzm1 ;
98 int nz_bc = nzm1 ;
99 if (par_bc != 0x0)
100 if (par_bc->get_n_int() > 0) nz_bc = par_bc->get_int() ;
101 n_shell = (nz_bc < n_shell ? nz_bc : n_shell) ;
102 bool cedbc = (ced && (nz_bc == nzm1)) ;
103#ifndef NDEBUG
104 if (!cedbc) {
105 assert(par_bc != 0x0) ;
106 assert(par_bc->get_n_tbl_mod() >= 3) ;
107 }
108#endif
109 int nt = mgrid.get_nt(0) ;
110 int np = mgrid.get_np(0) ;
111
112 Scalar source = aaa ;
113 Scalar source_coq = aaa ;
114 source_coq.annule_domain(0) ;
115 if (ced) source_coq.annule_domain(nzm1) ;
116 source_coq.mult_r() ;
117 source.set_spectral_va().ylm() ;
118 source_coq.set_spectral_va().ylm() ;
119 Base_val base = source.get_spectral_base() ;
120 base.mult_x() ;
121
122 tilde_vr.annule_hard() ;
123 tilde_vr.set_spectral_base(base) ;
124 tilde_vr.set_spectral_va().set_etat_cf_qcq() ;
125 tilde_vr.set_spectral_va().c_cf->annule_hard() ;
126 tilde_eta.annule_hard() ;
127 tilde_eta.set_spectral_base(base) ;
128 tilde_eta.set_spectral_va().set_etat_cf_qcq() ;
129 tilde_eta.set_spectral_va().c_cf->annule_hard() ;
130
131 Mtbl_cf sol_part_vr(mgrid, base) ; sol_part_vr.annule_hard() ;
132 Mtbl_cf sol_part_eta(mgrid, base) ; sol_part_eta.annule_hard() ;
133 Mtbl_cf sol_hom1_vr(mgrid, base) ; sol_hom1_vr.annule_hard() ;
134 Mtbl_cf sol_hom1_eta(mgrid, base) ; sol_hom1_eta.annule_hard() ;
135 Mtbl_cf sol_hom2_vr(mgrid, base) ; sol_hom2_vr.annule_hard() ;
136 Mtbl_cf sol_hom2_eta(mgrid, base) ; sol_hom2_eta.annule_hard() ;
137
138 int l_q, m_q, base_r ;
139
140 //---------------
141 //-- NUCLEUS ---
142 //---------------
143 {int lz = 0 ;
144 int nr = mgrid.get_nr(lz) ;
145 double alpha = mp_aff->get_alpha()[lz] ;
146 Matrice ope(2*nr, 2*nr) ;
147 ope.set_etat_qcq() ;
148
149 for (int k=0 ; k<np+1 ; k++) {
150 for (int j=0 ; j<nt ; j++) {
151 // quantic numbers and spectral bases
152 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
153 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 0)) {
154 Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
155 Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
156
157 for (int lin=0; lin<nr; lin++)
158 for (int col=0; col<nr; col++)
159 ope.set(lin,col) = md(lin,col) + 2*ms(lin,col) ;
160 for (int lin=0; lin<nr; lin++)
161 for (int col=0; col<nr; col++)
162 ope.set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
163 for (int lin=0; lin<nr; lin++)
164 for (int col=0; col<nr; col++)
165 ope.set(lin+nr,col) = -ms(lin,col) ;
166 for (int lin=0; lin<nr; lin++)
167 for (int col=0; col<nr; col++)
168 ope.set(lin+nr,col+nr) = md(lin,col) + ms(lin,col) ;
169
170 ope *= 1./alpha ;
171 int ind1 = nr ;
172
173 if (l_q==1) {
174 ind1 += 1 ;
175 int pari = 1 ;
176 for (int col=0 ; col<nr; col++) {
177 ope.set(nr-1,col) = pari ;
178 ope.set(nr-1,col+nr) = -pari ;
179 pari = - pari ;
180 }
181 ope.set(2*nr-1, nr) = 1;
182 }
183 else{
184 for (int col=0; col<2*nr; col++)
185 ope.set(ind1+nr-2, col) = 0 ;
186 for (int col=0; col<2*nr; col++) {
187 ope.set(nr-1, col) = 0 ;
188 ope.set(2*nr-1, col) = 0 ;
189 }
190 int pari = 1 ;
191 if (base_r == R_CHEBP) {
192 for (int col=0; col<nr; col++) {
193 ope.set(nr-1, col) = pari ;
194 ope.set(2*nr-1, col+nr) = pari ;
195 pari = - pari ;
196 }
197 }
198 else { //In the odd case, the last coefficient must be zero!
199 ope.set(nr-1, nr-1) = 1 ;
200 ope.set(2*nr-1, 2*nr-1) = 1 ;
201 }
202
203 ope.set(ind1+nr-2, ind1) = 1 ;
204 }
205
206 ope.set_lu() ;
207
208 Tbl sec(2*nr) ;
209 sec.set_etat_qcq() ;
210 for (int lin=0; lin<nr; lin++)
211 sec.set(lin) = 0 ;
212 for (int lin=0; lin<nr; lin++)
213 sec.set(nr+lin) = (*source.get_spectral_va().c_cf)
214 (lz, k, j, lin) ;
215 sec.set(2*nr-1) = 0 ;
216 sec.set(ind1+nr-2) = 0 ;
217 Tbl sol = ope.inverse(sec) ;
218
219
220
221 for (int i=0; i<nr; i++) {
222 sol_part_vr.set(lz, k, j, i) = sol(i) ;
223 sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
224 }
225 sec.annule_hard() ;
226 sec.set(ind1+nr-2) = 1 ;
227
228 sol = ope.inverse(sec) ;
229
230 for (int i=0; i<nr; i++) {
231 sol_hom2_vr.set(lz, k, j, i) = sol(i) ;
232 sol_hom2_eta.set(lz, k, j, i) = sol(i+nr) ;
233 }
234 }
235 }
236 }
237 }
238
239 //-------------
240 // -- Shells --
241 //-------------
242
243 for (int lz=1; lz <= n_shell; lz++) {
244 int nr = mgrid.get_nr(lz) ;
245 assert(mgrid.get_nt(lz) == nt) ;
246 assert(mgrid.get_np(lz) == np) ;
247 double alpha = mp_aff->get_alpha()[lz] ;
248 double ech = mp_aff->get_beta()[lz] / alpha ;
249 Matrice ope(2*nr, 2*nr) ;
250 ope.set_etat_qcq() ;
251
252 for (int k=0 ; k<np+1 ; k++) {
253 for (int j=0 ; j<nt ; j++) {
254 // quantic numbers and spectral bases
255 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
256 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 0)) {
257 Diff_xdsdx oxd(base_r, nr) ; const Matrice& mxd = oxd.get_matrice() ;
258 Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
259 Diff_id oid(base_r, nr) ; const Matrice& mid = oid.get_matrice() ;
260
261 for (int lin=0; lin<nr; lin++)
262 for (int col=0; col<nr; col++)
263 ope.set(lin,col) = mxd(lin,col) + ech*md(lin,col)
264 + 2*mid(lin,col) ;
265 for (int lin=0; lin<nr; lin++)
266 for (int col=0; col<nr; col++)
267 ope.set(lin,col+nr) = -l_q*(l_q+1)*mid(lin,col) ;
268 for (int lin=0; lin<nr; lin++)
269 for (int col=0; col<nr; col++)
270 ope.set(lin+nr,col) = -mid(lin,col) ;
271 for (int lin=0; lin<nr; lin++)
272 for (int col=0; col<nr; col++)
273 ope.set(lin+nr,col+nr) = mxd(lin,col) + ech*md(lin,col) + mid(lin,col) ;
274
275 int ind0 = 0 ;
276 int ind1 = nr ;
277 for (int col=0; col<2*nr; col++) {
278 ope.set(ind0+nr-1, col) = 0 ;
279 ope.set(ind1+nr-1, col) = 0 ;
280 }
281 ope.set(ind0+nr-1, ind0) = 1 ;
282 ope.set(ind1+nr-1, ind1) = 1 ;
283
284 ope.set_lu() ;
285
286 Tbl sec(2*nr) ;
287 sec.set_etat_qcq() ;
288 for (int lin=0; lin<nr; lin++)
289 sec.set(lin) = 0 ;
290 for (int lin=0; lin<nr; lin++)
291 sec.set(nr+lin) = (*source_coq.get_spectral_va().c_cf)
292 (lz, k, j, lin) ;
293 sec.set(ind0+nr-1) = 0 ;
294 sec.set(ind1+nr-1) = 0 ;
295 Tbl sol = ope.inverse(sec) ;
296 for (int i=0; i<nr; i++) {
297 sol_part_vr.set(lz, k, j, i) = sol(i) ;
298 sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
299 }
300
301
302 sec.annule_hard() ;
303 sec.set(ind0+nr-1) = 1 ;
304 sol = ope.inverse(sec) ;
305
306
307 for (int i=0; i<nr; i++) {
308 sol_hom1_vr.set(lz, k, j, i) = sol(i) ;
309 sol_hom1_eta.set(lz, k, j, i) = sol(i+nr) ;
310 }
311 sec.set(ind0+nr-1) = 0 ;
312 sec.set(ind1+nr-1) = 1 ;
313 sol = ope.inverse(sec) ;
314
315
316
317 for (int i=0; i<nr; i++) {
318 sol_hom2_vr.set(lz, k, j, i) = sol(i) ;
319 sol_hom2_eta.set(lz, k, j, i) = sol(i+nr) ;
320 }
321 }
322 }
323 }
324 }
325
326 //------------------------------
327 // Compactified external domain
328 //------------------------------
329 if (cedbc) {int lz = nzm1 ;
330 int nr = mgrid.get_nr(lz) ;
331 assert(mgrid.get_nt(lz) == nt) ;
332 assert(mgrid.get_np(lz) == np) ;
333 double alpha = mp_aff->get_alpha()[lz] ;
334 Matrice ope(2*nr, 2*nr) ;
335 ope.set_etat_qcq() ;
336
337 for (int k=0 ; k<np+1 ; k++) {
338 for (int j=0 ; j<nt ; j++) {
339 // quantic numbers and spectral bases
340 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
341 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 0)) {
342 Diff_dsdx od(base_r, nr) ; const Matrice& md = od.get_matrice() ;
343 Diff_sx os(base_r, nr) ; const Matrice& ms = os.get_matrice() ;
344
345
346 for (int lin=0; lin<nr; lin++)
347 for (int col=0; col<nr; col++)
348 ope.set(lin,col) = - md(lin,col) + 2*ms(lin,col) ;
349 for (int lin=0; lin<nr; lin++)
350 for (int col=0; col<nr; col++)
351 ope.set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
352 for (int lin=0; lin<nr; lin++)
353 for (int col=0; col<nr; col++)
354 ope.set(lin+nr,col) = -ms(lin,col) ;
355 for (int lin=0; lin<nr; lin++)
356 for (int col=0; col<nr; col++)
357 ope.set(lin+nr,col+nr) = -md(lin,col) + ms(lin,col) ;
358
359 ope *= 1./alpha ;
360 int ind0 = 0 ;
361 int ind1 = nr ;
362 for (int col=0; col<2*nr; col++) {
363 ope.set(ind0+nr-1, col) = 0 ;
364 ope.set(ind1+nr-2, col) = 0 ;
365 ope.set(ind1+nr-1, col) = 0 ;
366 }
367 for (int col=0; col<nr; col++) {
368 ope.set(ind0+nr-1, col+ind0) = 1 ;
369 ope.set(ind1+nr-1, col+ind1) = 1 ;
370 }
371 ope.set(ind1+nr-2, ind1+1) = 1 ;
372
373 ope.set_lu() ;
374
375 Tbl sec(2*nr) ;
376 sec.set_etat_qcq() ;
377 for (int lin=0; lin<nr; lin++)
378 sec.set(lin) = 0 ;
379 for (int lin=0; lin<nr; lin++)
380 sec.set(nr+lin) = (*source.get_spectral_va().c_cf)
381 (lz, k, j, lin) ;
382 sec.set(ind0+nr-1) = 0 ;
383 sec.set(ind1+nr-2) = 0 ;
384 sec.set(ind1+nr-1) = 0 ;
385 Tbl sol = ope.inverse(sec) ;
386
387 for (int i=0; i<nr; i++) {
388 sol_part_vr.set(lz, k, j, i) = sol(i) ;
389 sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
390 }
391 sec.annule_hard() ;
392 sec.set(ind1+nr-2) = 1 ;
393 sol = ope.inverse(sec) ;
394 for (int i=0; i<nr; i++) {
395 sol_hom1_vr.set(lz, k, j, i) = sol(i) ;
396 sol_hom1_eta.set(lz, k, j, i) = sol(i+nr) ;
397 }
398 }
399 }
400 }
401 }
402
403 int taille = 2*nz_bc + 1;
404 if (cedbc) taille-- ;
405 Mtbl_cf& mvr = *tilde_vr.set_spectral_va().c_cf ;
406 Mtbl_cf& meta = *tilde_eta.set_spectral_va().c_cf ;
407
408 Tbl sec_membre(taille) ;
409 Matrice systeme(taille, taille) ;
410 int ligne ; int colonne ;
411 Tbl pipo(1) ;
412 const Tbl& mub = (cedbc ? pipo : par_bc->get_tbl_mod(2) );
413 double c_vr = (cedbc ? 0 : par_bc->get_tbl_mod(0)(0) ) ;
414 double d_vr = (cedbc ? 0 : par_bc->get_tbl_mod(0)(1) ) ;
415 double c_eta = (cedbc ? 0 : par_bc->get_tbl_mod(0)(2) ) ;
416 double d_eta = (cedbc ? 0 : par_bc->get_tbl_mod(0)(3) ) ;
417
418
419
420 Mtbl_cf dhom1_vr = sol_hom1_vr ;
421 Mtbl_cf dhom2_vr = sol_hom2_vr ;
422 Mtbl_cf dpart_vr = sol_part_vr ;
423 Mtbl_cf dhom1_eta = sol_hom1_eta ;
424 Mtbl_cf dhom2_eta = sol_hom2_eta ;
425 Mtbl_cf dpart_eta = sol_part_eta ;
426 if (!cedbc) {
427 dhom1_vr.dsdx() ;
428 dhom2_vr.dsdx() ;
429 dpart_vr.dsdx() ;
430 dhom1_eta.dsdx() ;
431 dhom2_eta.dsdx() ;
432 dpart_eta.dsdx() ;
433 }
434
435 // Loop on l and m
436 //----------------
437 int nr ;
438 for (int k=0 ; k<np+1 ; k++)
439 for (int j=0 ; j<nt ; j++) {
440 base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
441 if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q > 0)) {
442 ligne = 0 ;
443 colonne = 0 ;
444 systeme.annule_hard() ;
445 sec_membre.annule_hard() ;
446
447
448 if ((nz==1)&&(mgrid.get_type_r(0) == RARE)) {
449 // Only one zone, which is a nucleus
450 double alpha = mp_aff->get_alpha()[nz_bc] ;
451 systeme.set(ligne, colonne) =
452 c_vr*sol_hom2_vr.val_out_bound_jk(nz_bc, j, k)
453 + d_vr*dhom2_vr.val_out_bound_jk(nz_bc, j, k) / alpha
454 + c_eta*sol_hom2_eta.val_out_bound_jk(nz_bc, j, k)
455 + d_eta*dhom2_eta.val_out_bound_jk(nz_bc, j, k) / alpha ;
456
457 sec_membre.set(ligne) -= c_vr*sol_part_vr.val_out_bound_jk(nz_bc, j, k)
458 + d_vr*dpart_vr.val_out_bound_jk(nz_bc, j, k)/alpha
459 + c_eta*sol_part_eta.val_out_bound_jk(nz_bc, j, k)
460 + d_eta*dpart_eta.val_out_bound_jk(nz_bc, j, k)/alpha
461 - mub(k, j) ;
462 }
463 else {
464 // General case, two zones at least
465 //Nucleus
466 systeme.set(ligne, colonne) = sol_hom2_vr.val_out_bound_jk(0, j, k) ;
467 sec_membre.set(ligne) = -sol_part_vr.val_out_bound_jk(0, j, k) ;
468 ligne++ ;
469
470 systeme.set(ligne, colonne) = sol_hom2_eta.val_out_bound_jk(0, j, k) ;
471 sec_membre.set(ligne) = -sol_part_eta.val_out_bound_jk(0, j, k) ;
472 colonne++ ;
473
474 //shells
475 for (int zone=1 ; zone<nz_bc ; zone++) {
476 nr = mgrid.get_nr(zone) ;
477 ligne-- ;
478
479 //Condition at x = -1
480 systeme.set(ligne, colonne) =
481 - sol_hom1_vr.val_in_bound_jk(zone, j, k) ;
482 systeme.set(ligne, colonne+1) =
483 - sol_hom2_vr.val_in_bound_jk(zone, j, k) ;
484
485 sec_membre.set(ligne) += sol_part_vr.val_in_bound_jk(zone, j, k) ;
486 ligne++ ;
487
488 systeme.set(ligne, colonne) =
489 - sol_hom1_eta.val_in_bound_jk(zone, j, k) ;
490 systeme.set(ligne, colonne+1) =
491 - sol_hom2_eta.val_in_bound_jk(zone, j, k) ;
492
493 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(zone, j, k) ;
494 ligne++ ;
495
496 // Condition at x=1
497 systeme.set(ligne, colonne) =
498 sol_hom1_vr.val_out_bound_jk(zone, j, k) ;
499 systeme.set(ligne, colonne+1) =
500 sol_hom2_vr.val_out_bound_jk(zone, j, k) ;
501
502 sec_membre.set(ligne) -= sol_part_vr.val_out_bound_jk(zone, j, k) ;
503 ligne++ ;
504
505 systeme.set(ligne, colonne) =
506 sol_hom1_eta.val_out_bound_jk(zone, j, k) ;
507 systeme.set(ligne, colonne+1) =
508 sol_hom2_eta.val_out_bound_jk(zone, j, k) ;
509
510 sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone, j, k) ;
511
512 colonne += 2 ;
513 }
514
515 //Last domain
516 nr = mgrid.get_nr(nz_bc) ;
517 double alpha = mp_aff->get_alpha()[nz_bc] ;
518 ligne-- ;
519
520 //Condition at x = -1
521 systeme.set(ligne, colonne) =
522 - sol_hom1_vr.val_in_bound_jk(nz_bc, j, k) ;
523 if (!cedbc) systeme.set(ligne, colonne+1) =
524 - sol_hom2_vr.val_in_bound_jk(nz_bc, j, k) ;
525
526 sec_membre.set(ligne) += sol_part_vr.val_in_bound_jk(nz_bc, j, k) ;
527 ligne++ ;
528
529 systeme.set(ligne, colonne) =
530 - sol_hom1_eta.val_in_bound_jk(nz_bc, j, k) ;
531 if (!cedbc) systeme.set(ligne, colonne+1) =
532 - sol_hom2_eta.val_in_bound_jk(nz_bc, j, k) ;
533
534 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(nz_bc, j, k) ;
535 ligne++ ;
536
537 if (!cedbc) {// Special condition at x=1
538 systeme.set(ligne, colonne) =
539 c_vr*sol_hom1_vr.val_out_bound_jk(nz_bc, j, k)
540 + d_vr*dhom1_vr.val_out_bound_jk(nz_bc, j, k) / alpha
541 + c_eta*sol_hom1_eta.val_out_bound_jk(nz_bc, j, k)
542 + d_eta*dhom1_eta.val_out_bound_jk(nz_bc, j, k) / alpha ;
543 systeme.set(ligne, colonne+1) =
544 c_vr*sol_hom2_vr.val_out_bound_jk(nz_bc, j, k)
545 + d_vr*dhom2_vr.val_out_bound_jk(nz_bc, j, k) / alpha
546 + c_eta*sol_hom2_eta.val_out_bound_jk(nz_bc, j, k)
547 + d_eta*dhom2_eta.val_out_bound_jk(nz_bc, j, k) / alpha ;
548
549 sec_membre.set(ligne) -= c_vr*sol_part_vr.val_out_bound_jk(nz_bc, j, k)
550 + d_vr*dpart_vr.val_out_bound_jk(nz_bc, j, k)/alpha
551 + c_eta*sol_part_eta.val_out_bound_jk(nz_bc, j, k)
552 + d_eta*dpart_eta.val_out_bound_jk(nz_bc, j, k)/alpha
553 - mub(k, j) ;
554 }
555 }
556
557 // Solution of the system giving the coefficients for the homogeneous
558 // solutions
559 //-------------------------------------------------------------------
560 systeme.set_lu() ;
561 Tbl facteur = systeme.inverse(sec_membre) ;
562 int conte = 0 ;
563
564
565 // everything is put to the right place...
566 //----------------------------------------
567 nr = mgrid.get_nr(0) ; //nucleus
568 for (int i=0 ; i<nr ; i++) {
569 mvr.set(0, k, j, i) = sol_part_vr(0, k, j, i)
570 + facteur(conte)*sol_hom2_vr(0, k, j, i) ;
571 meta.set(0, k, j, i) = sol_part_eta(0, k, j, i)
572 + facteur(conte)*sol_hom2_eta(0, k, j, i) ;
573 }
574 conte++ ;
575 for (int zone=1 ; zone<=n_shell ; zone++) { //shells
576 nr = mgrid.get_nr(zone) ;
577 for (int i=0 ; i<nr ; i++) {
578 mvr.set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
579 + facteur(conte)*sol_hom1_vr(zone, k, j, i)
580 + facteur(conte+1)*sol_hom2_vr(zone, k, j, i) ;
581
582 meta.set(zone, k, j, i) = sol_part_eta(zone, k, j, i)
583 + facteur(conte)*sol_hom1_eta(zone, k, j, i)
584 + facteur(conte+1)*sol_hom2_eta(zone, k, j, i) ;
585 }
586 conte+=2 ;
587 }
588 if (cedbc) {
589 nr = mgrid.get_nr(nzm1) ; //compactified external domain
590 for (int i=0 ; i<nr ; i++) {
591 mvr.set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
592 + facteur(conte)*sol_hom1_vr(nzm1, k, j, i) ;
593
594 meta.set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
595 + facteur(conte)*sol_hom1_eta(nzm1, k, j, i) ;
596 }
597 }
598
599 } // End of nullite_plm
600 } //End of loop on theta
601
602
603 if (tilde_vr.set_spectral_va().c != 0x0)
604 delete tilde_vr.set_spectral_va().c ;
605 tilde_vr.set_spectral_va().c = 0x0 ;
606 tilde_vr.set_spectral_va().ylm_i() ;
607
608 if (tilde_eta.set_spectral_va().c != 0x0)
609 delete tilde_eta.set_spectral_va().c ;
610 tilde_eta.set_spectral_va().c = 0x0 ;
611 tilde_eta.set_spectral_va().ylm_i() ;
612
613
614
615}
616}
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_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::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_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_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_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_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::get_n_int
int get_n_int() const
Returns the number of stored int 's addresses.
Definition param.C:242
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:393
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
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::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:375
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:364
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene::Valeur::set_etat_cf_qcq
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
Definition valeur.C:715
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:141
Lorene::Valeur::c
Mtbl * c
Values of the function at the points of the multi-grid.
Definition valeur.h:309
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:312
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm().
Definition valeur_ylm_i.C:134
Lorene::Vector_divfree::sol_Dirac_A
void sol_Dirac_A(const Scalar &aaa, Scalar &eta, Scalar &vr, const Param *par_bc=0x0) const
Solves a system of two-coupled first-order PDEs obtained from the divergence-free condition and the r...
Definition vector_divfree_A.C:80
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