LORENE
vector_df_poisson.C
1 /*
2  * Resolution of the divergence-free vector Poisson equation
3  *
4  * (see file vector.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2003 Eric Gourgoulhon & 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 vector_df_poisson_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/vector_df_poisson.C,v 1.15 2014/10/13 08:53:44 j_novak Exp $" ;
29 
30 /*
31  * $Id: vector_df_poisson.C,v 1.15 2014/10/13 08:53:44 j_novak Exp $
32  * $Log: vector_df_poisson.C,v $
33  * Revision 1.15 2014/10/13 08:53:44 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.14 2014/10/06 15:13:20 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.13 2005/02/15 09:45:22 j_novak
40  * Correction of an error in the matching.
41  *
42  * Revision 1.12 2005/02/09 16:53:11 j_novak
43  * Now V^r and eta are matched across domains, but not any of their derivatives.
44  *
45  * Revision 1.11 2005/02/09 14:52:01 j_novak
46  * Better solution in the shells.
47  *
48  * Revision 1.10 2005/02/09 13:20:27 j_novak
49  * Completely new way of solving the vector Poisson equation in the Div_free
50  * case: the system is inverted "as a whole" for V^r and eta. This works only
51  * with Map_af...
52  *
53  *
54  * $Header: /cvsroot/Lorene/C++/Source/Tensor/vector_df_poisson.C,v 1.15 2014/10/13 08:53:44 j_novak Exp $
55  *
56  */
57 
58 // C headers
59 #include <cassert>
60 #include <cmath>
61 
62 // Lorene headers
63 #include "tensor.h"
64 #include "diff.h"
65 #include "proto.h"
66 #include "param.h"
67 
68 namespace Lorene {
70 
71  // All this has a meaning only for spherical components:
72 #ifndef NDEBUG
73  const Base_vect_spher* bvs = dynamic_cast<const Base_vect_spher*>(triad) ;
74  assert(bvs != 0x0) ;
75 #endif
76 
77  int nitermax = par.get_int() ;
78  int& niter = par.get_int_mod() ;
79  double relax = par.get_double() ;
80  double precis = par.get_double(1) ;
81  Cmp& ss_khi = par.get_cmp_mod(0) ;
82  Cmp& ss_mu = par.get_cmp_mod(1) ;
83 
84  // Solution for the r-component
85  // ----------------------------
86 
87  Scalar source_r = *(cmp[0]) ;
88  source_r.mult_r() ;
89 
90  Param par_khi ;
91  par_khi.add_int(nitermax, 0) ;
92  par_khi.add_int_mod(niter, 0) ;
93  par_khi.add_double(relax, 0) ;
94  par_khi.add_double(precis, 1) ;
95  par_khi.add_cmp_mod(ss_khi, 0) ;
96 
97  Scalar khi (*mp) ;
98  khi.set_etat_zero() ;
99 
100  source_r.poisson(par_khi, khi) ;
101  khi.div_r() ; // khi now contains V^r
102 
103  // Solution for mu
104  // ---------------
105 
106  Param par_mu ;
107  par_mu.add_int(nitermax, 0) ;
108  par_mu.add_int_mod(niter, 0) ;
109  par_mu.add_double(relax, 0) ;
110  par_mu.add_double(precis, 1) ;
111  par_mu.add_cmp_mod(ss_mu, 0) ;
112 
113  Scalar mu_resu (*mp) ;
114  mu_resu.set_etat_zero() ;
115 
116  mu().poisson(par_mu, mu_resu) ;
117 
118  // Final result
119  // ------------
120 
121  Vector_divfree resu(*mp, *triad, *met_div) ;
122 
123  resu.set_vr_mu(khi, mu_resu) ;
124 
125  return resu ;
126 
127 }
128 
129 /*
130  * In the case without parameters, first is solved the equation for mu and then
131  * the system of equations for (eta, V^r) is inverted as a whole:
132  * d2 eta / dr2 + 2/r d eta / dr - 1/r d V^r / dr = S(eta)
133  * d V^r / dr + 2/r V^r - l(l+1)/r eta = 0 (div free condition)
134  *
135  * There is no l=0 contribution (divergence free in all space!).
136  * In the nucleus and the CED the system is inverted for h(r) and v(r) ,
137  * such that eta = r^2 h and V^r = r^2 v in the nucleus,
138  * in the compactified domain one has eta = u^2 h and V^r = u^2 v (where u=1/r);
139  * In the shells, both equations are multiplied by r.
140  * These methods are used only to get particular solutions.
141  *
142  * Homogeneous solutions are known analitically: r^(l-1) and/or 1/r^(l+2)
143  * It is then only possible to match eta and V^r, but not their derivatives,
144  * due to the div-free condition.
145  */
147 
148  const Map_af* mpaff = dynamic_cast<const Map_af*>(mp) ;
149 #ifndef NDEBUG
150  for (int i=0; i<3; i++)
151  assert(cmp[i]->check_dzpuis(4)) ;
152  // All this has a meaning only for spherical components:
153  const Base_vect_spher* bvs = dynamic_cast<const Base_vect_spher*>(triad) ;
154  assert(bvs != 0x0) ;
155  //## ... and affine mapping, for the moment!
156  assert( mpaff != 0x0) ;
157 #endif
158 
159  // Final result
160  // ------------
161  Vector_divfree resu(*mpaff, *triad, *met_div) ;
162 
163  // Solution for mu
164  // ---------------
165  Scalar mu_resu = mu().poisson() ;
166 
167  Scalar f_r(*mpaff) ;
168  if (cmp[0]->get_etat() == ETATZERO) { // no work needed ...
169  f_r.set_etat_zero() ;
170  resu.set_vr_mu(f_r, mu_resu) ;
171  return resu ;
172  }
173 
174  // Some working objects
175  //---------------------
176  const Mg3d& mg = *mpaff->get_mg() ;
177  int nz = mg.get_nzone() ; int nzm1 = nz - 1;
178  assert(mg.get_type_r(0) == RARE) ;
179  assert(mg.get_type_r(nzm1) == UNSURR) ;
180  Scalar S_r = *cmp[0] ;
181  Scalar S_eta = eta() ;
182  S_r.set_spectral_va().ylm() ;
183  S_eta.set_spectral_va().ylm() ;
184  const Base_val& base = S_eta.get_spectral_va().base ;
185  Mtbl_cf sol_part_eta(mg, base) ; sol_part_eta.annule_hard() ;
186  Mtbl_cf sol_part_vr(mg, base) ; sol_part_vr.annule_hard() ;
187  Mtbl_cf solution_hom_un(mg, base) ; solution_hom_un.annule_hard() ;
188  Mtbl_cf solution_hom_deux(mg, base) ; solution_hom_deux.annule_hard() ;
189 
190  // Build-up & inversion of the system for (eta, V^r) in each domain
191  //-----------------------------------------------------------------
192 
193  // Nucleus
194  //--------
195  int nr = mg.get_nr(0) ;
196  int nt = mg.get_nt(0) ;
197  int np = mg.get_np(0) ;
198  double alpha = mpaff->get_alpha()[0] ;
199  double beta = mpaff->get_beta()[0] ;
200  int l_q = 0 ; int m_q = 0; int base_r = 0 ;
201  int nr0 = nr - 1 ; //one degree of freedom less because of division by r^2
202 
203  // Loop on l and m
204  //----------------
205  for (int k=0 ; k<np+1 ; k++) {
206  for (int j=0 ; j<nt ; j++) {
207  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
208  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
209  int dege1 = (l_q <3 ? 0 : 1) ; //degeneracy of eq.1
210  int dege2 = 0 ; //degeneracy of eq.2
211  int nr_eq1 = nr0 - dege1 ; //Eq.1 is for h (eta/r^2)
212  int nr_eq2 = nr0 - dege2 ; //Eq.2 is the div-free condition
213  int nrtot = nr_eq1 + nr_eq2 ;
214  Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
215  Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
216  Diff_x2dsdx2 d2(base_r, nr) ; const Matrice& md2 = d2.get_matrice() ;
217  Diff_xdsdx xd(base_r, nr) ; const Matrice& mxd = xd.get_matrice() ;
218  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
219 
220  // Building the operator
221  //----------------------
222  for (int lin=0; lin<nr_eq1; lin++) { //eq.1
223  for (int col=dege1; col<nr0; col++)
224  oper.set(lin,col-dege1)
225  = md2(lin,col) + 6*mxd(lin,col) + 6*mid(lin,col) ;
226  for (int col=dege2; col<nr0; col++)
227  oper.set(lin,col-dege2+nr_eq1) = -mxd(lin,col) -2*mid(lin,col) ;
228  }
229  for (int lin=0; lin<nr0; lin++) { //eq.2
230  for (int col=dege1; col<nr0; col++)
231  oper.set(lin+nr_eq1,col-dege1) = -l_q*(l_q+1)*mid(lin,col) ;
232  for (int col=dege2; col<nr0; col++)
233  oper.set(lin+nr_eq1, col-dege2+nr_eq1) = mxd(lin,col) + 4*mid(lin,col) ;
234  }
235  oper.set_lu() ;
236 
237  // Filling the r.h.s
238  //------------------
239  for (int i=0; i<nr_eq1; i++) //eq.1
240  sec_membre.set(i) = (*S_eta.get_spectral_va().c_cf)(0, k, j, i) ;
241  for (int i=0; i<nr0; i++) //eq.2
242  sec_membre.set(i+nr_eq1) = 0 ;
243 
244  // Inversion of the "big" operator
245  //--------------------------------
246  Tbl big_res = oper.inverse(sec_membre) ;
247  Tbl res_eta(nr) ; res_eta.set_etat_qcq() ;
248  Tbl res_vr(nr) ; res_vr.set_etat_qcq() ;
249 
250  // Putting coefficients of h and v to individual arrays
251  //-----------------------------------------------------
252  for (int i=0; i<dege1; i++)
253  res_eta.set(i) = 0 ;
254  for (int i=dege1; i<nr0; i++)
255  res_eta.set(i) = big_res(i-dege1) ;
256  res_eta.set(nr0) = 0 ;
257  for (int i=0; i<dege2; i++)
258  res_vr.set(i) = 0 ;
259  for (int i=dege2; i<nr0; i++)
260  res_vr.set(i) = big_res(i-dege2+nr_eq1) ;
261  res_vr.set(nr0) = 0 ;
262 
263  // Multiplication by xi^2 (the alpha^2 factor comes soon)
264  //-------------------------------------------------------
265  multx2_1d(nr, &res_eta.t, base_r) ;
266  multx2_1d(nr, &res_vr.t, base_r) ;
267 
268  // Homogeneous solution (only r^(l-1) in the nucleus)
269  Tbl sol_hom = solh(nr, l_q-1, 0., base_r) ;
270  for (int i=0 ; i<nr ; i++) {
271  sol_part_eta.set(0, k, j, i) = alpha*alpha*res_eta(i) ;
272  sol_part_vr.set(0, k, j, i) = alpha*alpha*res_vr(i) ;
273  solution_hom_un.set(0, k, j, i) = sol_hom(i) ;
274  solution_hom_deux.set(0, k, j, i) = 0. ;
275  }
276  }
277  }
278  }
279 
280 
281  // Shells
282  //-------
283  for (int zone=1 ; zone<nzm1 ; zone++) {
284  nr = mg.get_nr(zone) ;
285  assert (nr > 5) ;
286  assert(nt == mg.get_nt(zone)) ;
287  assert(np == mg.get_np(zone)) ;
288  alpha = mpaff->get_alpha()[zone] ;
289  beta = mpaff->get_beta()[zone] ;
290  double ech = beta / alpha ;
291 
292  // Loop on l and m
293  //----------------
294  for (int k=0 ; k<np+1 ; k++) {
295  for (int j=0 ; j<nt ; j++) {
296  base.give_quant_numbers(zone, k, j, m_q, l_q, base_r) ;
297  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
298  int dege1 = 2 ; //degeneracy of eq.1
299  int dege2 = 0 ; //degeneracy of eq.2
300  int nr_eq1 = nr - dege1 ; //Eq.1 is for eta
301  int nr_eq2 = nr - dege2 ; //Eq.2 is the div-free condition
302  int nrtot = nr_eq1 + nr_eq2 + 1;
303  Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
304  Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
305  Diff_x2dsdx2 x2d2(base_r, nr+1); const Matrice& m2d2 = x2d2.get_matrice() ;
306  Diff_xdsdx2 xd2(base_r, nr+1) ; const Matrice& mxd2 = xd2.get_matrice() ;
307  Diff_dsdx2 d2(base_r, nr+1) ; const Matrice& md2 = d2.get_matrice() ;
308  Diff_xdsdx xd(base_r, nr+1) ; const Matrice& mxd = xd.get_matrice() ;
309  Diff_dsdx d1(base_r, nr+1) ; const Matrice& md = d1.get_matrice() ;
310  Diff_id id(base_r, nr+1) ; const Matrice& mid = id.get_matrice() ;
311 
312  // Building the operator
313  //----------------------
314  for (int lin=0; lin<nr_eq1; lin++) {
315  for (int col=dege1; col<nr; col++)
316  oper.set(lin,col-dege1)
317  = mxd2(lin,col) + ech*md2(lin,col) + 2*md(lin,col) ;
318  for (int col=dege2; col<nr+1; col++)
319  oper.set(lin,col-dege2+nr_eq1) = -md(lin,col) ;
320  }
321  for (int lin=0; lin<nr_eq2; lin++) {
322  for (int col=dege1; col<nr; col++)
323  oper.set(lin+nr_eq1,col-dege1) = -l_q*(l_q+1)*mid(lin,col) ;
324  for (int col=dege2; col<nr+1; col++)
325  oper.set(lin+nr_eq1, col-dege2+nr_eq1) =
326  mxd(lin,col) + ech*md(lin,col) + 2*mid(lin,col) ;
327  }
328  //Additional line to avoid spurious homogeneous solutions
329  //this line is the first one of the V^r eq.
330  for (int col=dege1; col<nr; col++)
331  oper.set(nrtot-1, col-dege1) = 0 ;
332  for (int col=dege2; col<nr+1; col++)
333  oper.set(nrtot-1, col-dege2+nr_eq1) =
334  m2d2(0,col) + ech*(2*mxd2(0,col) + ech*md2(0,col))
335  +4*(mxd(0,col) +ech*md(0,col))
336  +(2 - l_q*(l_q+1))*mid(0,col) ;
337  oper.set_lu() ;
338 
339  // Filling the r.h.s
340  //------------------
341  Tbl sr(5) ; sr.set_etat_qcq() ;
342  Tbl seta(nr) ; seta.set_etat_qcq() ;
343  for (int i=0; i<5; i++) {
344  sr.set(i) = (*S_r.get_spectral_va().c_cf)(zone, k, j, i);
345  seta.set(i) = (*S_eta.get_spectral_va().c_cf)(zone, k, j, i) ;
346  }
347  for (int i=5; i<nr; i++)
348  seta.set(i) = (*S_eta.get_spectral_va().c_cf)(zone, k, j, i) ;
349  Tbl xsr= sr ; Tbl x2sr= sr ;
350  Tbl xseta= seta ;
351  multx2_1d(5, &x2sr.t, base_r) ; multx_1d(5, &xsr.t, base_r) ;
352  multx_1d(nr, &xseta.t, base_r) ;
353 
354  for (int i=0; i<nr_eq1; i++)
355  sec_membre.set(i) = alpha*(alpha*xseta(i) + beta*seta(i)) ;
356  for (int i=0; i<nr_eq2; i++)
357  sec_membre.set(i+nr_eq1) = 0 ;
358  sec_membre.set(nr_eq1+nr_eq2) = alpha*alpha*x2sr(0) + 2*alpha*beta*xsr(0)
359  + beta*beta*sr(0) ;
360 
361  // Inversion of the "big" operator
362  //--------------------------------
363  Tbl big_res = oper.inverse(sec_membre) ;
364  Tbl res_eta(nr) ; res_eta.set_etat_qcq() ;
365  Tbl res_vr(nr) ; res_vr.set_etat_qcq() ;
366 
367  // Putting coefficients of h and v to individual arrays
368  //-----------------------------------------------------
369  for (int i=0; i<dege1; i++)
370  res_eta.set(i) = 0 ;
371  for (int i=dege1; i<nr; i++)
372  res_eta.set(i) = big_res(i-dege1) ;
373  for (int i=0; i<dege2; i++)
374  res_vr.set(i) = 0 ;
375  for (int i=dege2; i<nr; i++)
376  res_vr.set(i) = big_res(i-dege2+nr_eq1) ;
377 
378  //homogeneous solutions
379  Tbl sol_hom1 = solh(nr, l_q-1, ech, base_r) ;
380  Tbl sol_hom2 = solh(nr, l_q+1, ech, base_r) ;
381  for (int i=0 ; i<nr ; i++) {
382  sol_part_eta.set(zone, k, j, i) = res_eta(i) ;
383  sol_part_vr.set(zone, k, j, i) = res_vr(i) ;
384  solution_hom_un.set(zone, k, j, i) = sol_hom1(0,i) ;
385  solution_hom_deux.set(zone, k, j, i) = sol_hom2(1,i) ;
386  }
387  }
388  }
389  }
390  }
391 
392  // Compactified external domain
393  //-----------------------------
394  nr = mg.get_nr(nzm1) ;
395  assert(nt == mg.get_nt(nzm1)) ;
396  assert(np == mg.get_np(nzm1)) ;
397  alpha = mpaff->get_alpha()[nzm1] ;
398  assert (nr > 4) ;
399  nr0 = nr - 2 ; //two degrees of freedom less because of division by r^2
400 
401  // Loop on l and m
402  //----------------
403  for (int k=0 ; k<np+1 ; k++) {
404  for (int j=0 ; j<nt ; j++) {
405  base.give_quant_numbers(nzm1, k, j, m_q, l_q, base_r) ;
406  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
407  int dege1 = 0; //degeneracy of eq.1
408  int dege2 = 1; //degeneracy of eq.2
409  int nr_eq1 = nr0 - dege1 ; //Eq.1 is for eta
410  int nr_eq2 = nr0 - dege2 ; //Eq.2 is the div-free condition
411  int nrtot = nr_eq1 + nr_eq2 ;
412  Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
413  Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
414  Diff_x2dsdx2 x2d2(base_r, nr) ; const Matrice& m2d2 = x2d2.get_matrice() ;
415  Diff_xdsdx xd(base_r, nr) ; const Matrice& mxd = xd.get_matrice() ;
416  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
417 
418  // Building the operator
419  //----------------------
420  for (int lin=0; lin<nr_eq1; lin++) {
421  for (int col=dege1; col<nr0; col++)
422  oper.set(lin,col-dege1)
423  = m2d2(lin,col) + 4*mxd(lin,col) + 2*mid(lin,col) ;
424  for (int col=dege2; col<nr0; col++)
425  oper.set(lin,col-dege2+nr_eq1) =
426  mxd(lin,col) + 2*mid(lin,col) ;
427  }
428  for (int lin=0; lin<nr_eq2; lin++) {
429  for (int col=dege1; col<nr0; col++)
430  oper.set(lin+nr_eq1,col-dege1) = l_q*(l_q+1)*mid(lin,col) ;
431  for (int col=dege2; col<nr0; col++)
432  oper.set(lin+nr_eq1, col-dege2+nr_eq1) = mxd(lin,col) ;
433  }
434  oper.set_lu() ;
435 
436  // Filling the r.h.s
437  //------------------
438  for (int i=0; i<nr_eq1; i++)
439  sec_membre.set(i) = (*S_eta.get_spectral_va().c_cf)(nzm1, k, j, i) ;
440  for (int i=0; i<nr_eq2; i++)
441  sec_membre.set(i+nr_eq1) = 0 ;
442  Tbl big_res = oper.inverse(sec_membre) ;
443  Tbl res_eta(nr) ; res_eta.set_etat_qcq() ;
444  Tbl res_vr(nr) ; res_vr.set_etat_qcq() ;
445 
446  // Putting coefficients of h and v to individual arrays
447  //-----------------------------------------------------
448  for (int i=0; i<dege1; i++)
449  res_eta.set(i) = 0 ;
450  for (int i=dege1; i<nr0; i++)
451  res_eta.set(i) = big_res(i-dege1) ;
452  res_eta.set(nr0) = 0 ;
453  res_eta.set(nr0+1) = 0 ;
454  for (int i=0; i<dege2; i++)
455  res_vr.set(i) = 0 ;
456  for (int i=dege2; i<nr0; i++)
457  res_vr.set(i) = big_res(i-dege2+nr_eq1) ;
458  res_vr.set(nr0) = 0 ;
459  res_vr.set(nr0+1) = 0 ;
460 
461  // Multiplication by r^2
462  //-----------------------
463  Tbl res_v2(nr) ; res_v2.set_etat_qcq() ;
464  Tbl res_e2(nr) ; res_e2.set_etat_qcq() ;
465  mult2_xm1_1d_cheb(nr, res_eta.t, res_e2.t) ;
466  mult2_xm1_1d_cheb(nr, res_vr.t, res_v2.t) ;
467 
468  // Homogeneous solution (only 1/r^(l+2) in the CED)
469  Tbl sol_hom = solh(nr, l_q+1, 0., base_r) ;
470  for (int i=0 ; i<nr ; i++) {
471  sol_part_eta.set(nzm1, k, j, i) = alpha*alpha*res_e2(i) ;
472  sol_part_vr.set(nzm1, k, j, i) = alpha*alpha*res_v2(i) ;
473  solution_hom_un.set(nzm1, k, j, i) = sol_hom(i) ;
474  solution_hom_deux.set(nzm1, k, j, i) = 0. ;
475  }
476  }
477  }
478  }
479 
480  // Now let's match everything ...
481  //-------------------------------
482 
483  // Resulting V^r & eta
484  Scalar vr(*mpaff) ; vr.set_etat_qcq() ;
485  vr.set_spectral_base(base) ;
487  Mtbl_cf& cf_vr = *vr.set_spectral_va().c_cf ;
488  cf_vr.annule_hard() ;
489  Scalar het(*mpaff) ; het.set_etat_qcq() ;
490  het.set_spectral_base(base) ;
492  Mtbl_cf& cf_eta = *het.set_spectral_va().c_cf ;
493  cf_eta.annule_hard() ;
494  int taille = 2*nzm1 ;
495  Tbl sec_membre(taille) ;
496  Matrice systeme(taille, taille) ;
497  systeme.set_etat_qcq() ;
498  int ligne ; int colonne ;
499 
500  // Loop on l and m
501  //----------------
502  for (int k=0 ; k<np+1 ; k++)
503  for (int j=0 ; j<nt ; j++) {
504  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
505  if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q != 0)) {
506 
507  ligne = 0 ;
508  colonne = 0 ;
509  sec_membre.annule_hard() ;
510  for (int l=0; l<taille; l++)
511  for (int c=0; c<taille; c++)
512  systeme.set(l,c) = 0 ;
513  //Nucleus
514  nr = mg.get_nr(0) ;
515  alpha = mpaff->get_alpha()[0] ;
516  // value of x^(l-1) at 1 ...
517  systeme.set(ligne, colonne) = 1. ;
518  for (int i=0 ; i<nr ; i++)
519  sec_membre.set(ligne) -= sol_part_eta(0, k, j, i) ;
520  ligne++ ;
521  // ... and of its couterpart for V^r
522  systeme.set(ligne, colonne) = l_q;
523  for (int i=0; i<nr; i++)
524  sec_membre.set(ligne) -= sol_part_vr(0,k,j,i) ;
525  colonne++ ;
526  //shells
527  for (int zone=1 ; zone<nzm1 ; zone++) {
528  nr = mg.get_nr(zone) ;
529  alpha = mpaff->get_alpha()[zone] ;
530  double echelle = mpaff->get_beta()[zone]/alpha ;
531  ligne -- ;
532  //value of (x+echelle)^(l-1) at -1
533  systeme.set(ligne, colonne) = -pow(echelle-1., double(l_q-1)) ;
534  // value of 1/(x+echelle) ^(l+2) at -1
535  systeme.set(ligne, colonne+1) = -1/pow(echelle-1., double(l_q+2)) ;
536  for (int i=0 ; i<nr ; i++)
537  if (i%2 == 0)
538  sec_membre.set(ligne) += sol_part_eta(zone, k, j, i) ;
539  else sec_membre.set(ligne) -= sol_part_eta(zone, k, j, i) ;
540  ligne++ ;
541  // ... and their couterparts for V^r
542  systeme.set(ligne, colonne) = -l_q*pow(echelle-1., double(l_q-1)) ;
543  systeme.set(ligne, colonne+1) = (l_q+1)/pow(echelle-1., double(l_q+2));
544  for (int i=0 ; i<nr ; i++)
545  if (i%2 == 0)
546  sec_membre.set(ligne) += sol_part_vr(zone, k, j, i) ;
547  else sec_membre.set(ligne) -= sol_part_vr(zone, k, j, i) ;
548  ligne++ ;
549 
550  //value of (x+echelle)^(l-1) at 1 :
551  systeme.set(ligne, colonne) = pow(echelle+1., double(l_q-1)) ;
552  // value of 1/(x+echelle)^(l+2) at 1 :
553  systeme.set(ligne, colonne+1) = 1./pow(echelle+1., double(l_q+2)) ;
554  for (int i=0 ; i<nr ; i++)
555  sec_membre.set(ligne) -= sol_part_eta(zone, k, j, i) ;
556  ligne ++ ;
557  //... and their couterparts for V^r
558  systeme.set(ligne, colonne) = l_q*pow(echelle+1., double(l_q-1)) ;
559  systeme.set(ligne, colonne+1) = -double(l_q+1)
560  / pow(echelle+1., double(l_q+2)) ;
561  for (int i=0 ; i<nr ; i++)
562  sec_membre.set(ligne) -= sol_part_vr(zone, k, j, i);
563  colonne += 2 ;
564  }
565  //Compactified external domain
566  nr = mg.get_nr(nzm1) ;
567 
568  alpha = mpaff->get_alpha()[nzm1] ;
569  ligne -- ;
570  //value of (x-1)^(l+2) at -1 :
571  systeme.set(ligne, colonne) = -pow(-2, double(l_q+2)) ;
572  for (int i=0 ; i<nr ; i++)
573  if (i%2 == 0) sec_membre.set(ligne) += sol_part_eta(nzm1, k, j, i) ;
574  else sec_membre.set(ligne) -= sol_part_eta(nzm1, k, j, i) ;
575  //... and of its couterpart for V^r
576  systeme.set(ligne+1, colonne) = double(l_q+1)*pow(-2, double(l_q+2)) ;
577  for (int i=0 ; i<nr ; i++)
578  if (i%2 == 0) sec_membre.set(ligne+1) += sol_part_vr(nzm1, k, j, i) ;
579  else sec_membre.set(ligne+1) -= sol_part_vr(nzm1, k, j, i) ;
580 
581  // Solution of the system giving the coefficients for the homogeneous
582  // solutions
583  //-------------------------------------------------------------------
584  if (taille > 2) systeme.set_band(2,2) ;
585  else systeme.set_band(1,1) ;
586  systeme.set_lu() ;
587  Tbl facteurs(systeme.inverse(sec_membre)) ;
588  int conte = 0 ;
589 
590  // everything is put to the right place, the same combination of hom.
591  // solutions (with some l or -(l+1) factors) must be used for V^r
592  //-------------------------------------------------------------------
593  nr = mg.get_nr(0) ; //nucleus
594  for (int i=0 ; i<nr ; i++) {
595  cf_eta.set(0, k, j, i) = sol_part_eta(0, k, j, i)
596  +facteurs(conte)*solution_hom_un(0, k, j, i) ;
597  cf_vr.set(0, k, j, i) = sol_part_vr(0, k, j, i)
598  +double(l_q)*facteurs(conte)*solution_hom_un(0, k, j, i) ;
599  }
600  conte++ ;
601  for (int zone=1 ; zone<nzm1 ; zone++) { //shells
602  nr = mg.get_nr(zone) ;
603  for (int i=0 ; i<nr ; i++) {
604  cf_eta.set(zone, k, j, i) =
605  sol_part_eta(zone, k, j, i)
606  +facteurs(conte)*solution_hom_un(zone, k, j, i)
607  +facteurs(conte+1)*solution_hom_deux(zone, k, j, i) ;
608  cf_vr.set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
609  +double(l_q)*facteurs(conte)*solution_hom_un(zone, k, j, i)
610  -double(l_q+1)*facteurs(conte+1)*solution_hom_deux(zone, k, j, i) ;
611  }
612  conte+=2 ;
613  }
614  nr = mg.get_nr(nz-1) ; //compactified external domain
615  for (int i=0 ; i<nr ; i++) {
616  cf_eta.set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
617  +facteurs(conte)*solution_hom_un(nzm1, k, j, i) ;
618  cf_vr.set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
619  -double(l_q+1)*facteurs(conte)*solution_hom_un(nzm1, k, j, i) ;
620  }
621  } // End of nullite_plm
622  } //End of loop on theta
623 
624  vr.set_spectral_va().ylm_i() ;
625  het.set_spectral_va().ylm_i() ;
626 
627  resu.set_vr_eta_mu(vr, het, mu_resu) ;
628 
629  return resu ;
630 
631 }
632 
633 }
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_x2dsdx2.C:95
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:443
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:172
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:477
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:246
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
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
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:324
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition: matrice.C:392
Lorene prototypes.
Definition: app_hor.h:64
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:294
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition: matrice.C:424
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_dsdx.C:94
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: scalar_pde.C:136
const Metric *const met_div
Metric with respect to which the divergence is defined.
Definition: vector.h:730
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tensor.h:303
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:129
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:353
Cmp & get_cmp_mod(int position=0) const
Returns the reference of a modifiable Cmp stored in the list.
Definition: param.C:1049
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Class for the elementary differential operator Identity (see the base class Diff ).
Definition: diff.h:210
void div_r()
Division by r everywhere; dzpuis is not changed.
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_xdsdx2.C:97
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:490
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:292
Matrix handling.
Definition: matrice.h:152
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition: tensor.h:315
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_dsdx2.C:91
double * t
The array of double.
Definition: tbl.h:173
const double * get_beta() const
Returns the pointer on the array beta.
Definition: map_af.C:481
Parameter storage.
Definition: param.h:125
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Definition: param.C:430
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Vector_divfree poisson() const
Computes the solution of a vectorial Poisson equation with *this as a source: .
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written: .
double & set(int j, int i)
Read/write of a particuliar element.
Definition: matrice.h:277
void set_band(int up, int low) const
Calculate the band storage of *std.
Definition: matrice.C:364
Spherical orthonormal vectorial bases (triads).
Definition: base_vect.h:308
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va.
Definition: scalar.C:797
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Multi-domain grid.
Definition: grilles.h:273
Bases of the spectral expansions.
Definition: base_val.h:322
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition: mtbl_cf.C:312
Affine radial mapping.
Definition: map.h:2027
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:409
void set_vr_mu(const Scalar &vr_i, const Scalar &mu_i)
Sets the angular potentials (see member p_mu ), and the component of the vector.
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:315
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:186
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written: .
Definition: vector_etamu.C:98
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:175
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:531
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition: param.C:361
Basic array class.
Definition: tbl.h:161
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
Definition: param.C:1004
Divergence-free vectors.
Definition: vector.h:724
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:474
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:295
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:372
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_xdsdx.C:98
void set_vr_eta_mu(const Scalar &vr_i, const Scalar &eta_i, const Scalar &mu_i)
Defines the components through , and .
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition: param.C:385
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:601