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