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 char 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 
76 namespace Lorene {
77 void Vector_divfree::sol_Dirac_A(const Scalar& aaa, Scalar& tilde_vr, Scalar& tilde_eta,
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() ;
98  n_shell = (nz_bc < n_shell ? nz_bc : n_shell) ;
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 ;
110  Scalar source_coq = 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 
417  Mtbl_cf dhom1_vr = sol_hom1_vr ;
418  Mtbl_cf dhom2_vr = sol_hom2_vr ;
419  Mtbl_cf dpart_vr = sol_part_vr ;
420  Mtbl_cf dhom1_eta = sol_hom1_eta ;
421  Mtbl_cf dhom2_eta = sol_hom2_eta ;
422  Mtbl_cf dpart_eta = sol_part_eta ;
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)
576  + facteur(conte)*sol_hom1_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)
580  + facteur(conte)*sol_hom1_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)
592  + facteur(conte)*sol_hom1_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 }
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Lorene::Base_val
Bases of the spectral expansions.
Definition: base_val.h:322
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:462
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:712
Lorene::Matrice::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:175
Lorene::Scalar::set_spectral_base
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
Definition: scalar.C:797
Lorene::Mg3d
Multi-domain grid.
Definition: grilles.h:273
Lorene::Diff_sx
Class for the elementary differential operator division by (see the base class Diff ).
Definition: diff.h:329
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
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:457
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:372
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
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:474
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:636
Lorene::Diff_id::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_id.C:91
Lorene::Matrice::annule_hard
void annule_hard()
Sets the logical state to ETATQCQ (undefined state).
Definition: matrice.C:193
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:666
Lorene::Mtbl_cf::set
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:294
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Lorene::Valeur::c
Mtbl * c
Values of the function at the points of the multi-grid
Definition: valeur.h:299
Lorene::Scalar::get_spectral_base
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va
Definition: scalar.h:1294
Lorene::Diff_sx::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_sx.C:100
Lorene::Scalar::get_spectral_va
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:601
Lorene::Matrice::inverse
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition: matrice.C:424
Lorene::Param::get_n_int
int get_n_int() const
Returns the number of stored int 's addresses.
Definition: param.C:239
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:452
Lorene::Diff_xdsdx::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_xdsdx.C:98
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
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:81
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Lorene::Scalar::mult_r
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
Definition: scalar_r_manip.C:208
Lorene::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:295
Lorene::Mtbl_cf::annule_hard
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition: mtbl_cf.C:312
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:77
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Diff_xdsdx
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:409
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:452
Lorene::Map_af
Affine radial mapping.
Definition: map.h:2027
Lorene::Diff_id
Class for the elementary differential operator Identity (see the base class Diff ).
Definition: diff.h:210
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:428
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
Lorene::Matrice
Matrix handling.
Definition: matrice.h:152
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:584
Lorene::Map_af::get_beta
const double * get_beta() const
Returns the pointer on the array beta.
Definition: map_af.C:481
Lorene::Param
Parameter storage.
Definition: param.h:125
Lorene::Base_val::mult_x
void mult_x()
The basis is transformed as with a multiplication by .
Definition: base_val_manip.C:157
Lorene::Scalar::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:380
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:292
Lorene::Diff_dsdx::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_dsdx.C:94
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
Lorene::Matrice::set
double & set(int j, int i)
Read/write of a particuliar element.
Definition: matrice.h:277
Lorene::Diff_dsdx
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:129
Lorene::Map_af::get_alpha
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:477
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:186
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131
Lorene::Matrice::set_lu
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition: matrice.C:392
Lorene::Mtbl_cf::dsdx
void dsdx()
Definition: valeur_dsdx.C:153