LORENE
vector_poisson_boundary.C
1 /*
2  * Method for vector Poisson equation inverting eqs. for V^r and eta as a block
3  * (with a boundary condition on the excised sphere).
4  *
5  * (see file vector.h for documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 2005 Jerome Novak
11  * Francois Limousin
12  * Jose Luis Jaramillo
13  *
14  * This file is part of LORENE.
15  *
16  * LORENE is free software; you can redistribute it and/or modify
17  * it under the terms of the GNU General Public License version 2
18  * as published by the Free Software Foundation.
19  *
20  * LORENE is distributed in the hope that it will be useful,
21  * but WITHOUT ANY WARRANTY; without even the implied warranty of
22  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23  * GNU General Public License for more details.
24  *
25  * You should have received a copy of the GNU General Public License
26  * along with LORENE; if not, write to the Free Software
27  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28  *
29  */
30 
31 
32 
33 /*
34  * $Id: vector_poisson_boundary.C,v 1.4 2016/12/05 16:18:18 j_novak Exp $
35  * $Log: vector_poisson_boundary.C,v $
36  * Revision 1.4 2016/12/05 16:18:18 j_novak
37  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
38  *
39  * Revision 1.3 2014/10/13 08:53:45 j_novak
40  * Lorene classes and functions now belong to the namespace Lorene.
41  *
42  * Revision 1.2 2014/10/06 15:13:21 j_novak
43  * Modified #include directives to use c++ syntax.
44  *
45  * Revision 1.1 2005/06/09 08:00:09 f_limousin
46  * Implement a new function sol_elliptic_boundary() and
47  * Vector::poisson_boundary(...) which solve the vectorial poisson
48  * equation (method 6) with an inner boundary condition.
49  *
50  *
51  * $Header: /cvsroot/Lorene/C++/Source/Tensor/vector_poisson_boundary.C,v 1.4 2016/12/05 16:18:18 j_novak Exp $
52  *
53  */
54 
55 // C headers
56 #include <cassert>
57 #include <cstdlib>
58 #include <cmath>
59 
60 // Lorene headers
61 #include "metric.h"
62 #include "diff.h"
63 #include "param_elliptic.h"
64 #include "proto.h"
65 #include "utilitaires.h"
66 
67 namespace Lorene {
68 void Vector::poisson_boundary(double lam, const Mtbl_cf& bound_vr,
69  const Mtbl_cf& bound_eta, const Mtbl_cf& bound_mu,
70  int num_front, double fact_dir, double fact_neu,
71  Vector& resu) const {
72 
73  const Map_af* mpaff = dynamic_cast<const Map_af*>(mp) ;
74 #ifndef NDEBUG
75  for (int i=0; i<3; i++)
76  assert(cmp[i]->check_dzpuis(4)) ;
77  // All this has a meaning only for spherical components:
78  const Base_vect_spher* bvs = dynamic_cast<const Base_vect_spher*>(triad) ;
79  assert(bvs != 0x0) ;
80  //## ... and affine mapping, for the moment!
81  assert( mpaff != 0x0) ;
82 #endif
83 
84  if (fabs(lam + 1.) < 1.e-8) {
85  cout << "Not implemented yet !!" << endl ;
86  abort() ;
87  /*
88  const Metric_flat& mets = mp->flat_met_spher() ;
89  Vector_divfree sou(*mp, *triad, mets) ;
90  for (int i=1; i<=3; i++) sou.set(i) = *cmp[i-1] ;
91  resu = sou.poisson() ;
92  return ;
93  */
94  }
95 
96  // Some working objects
97  //---------------------
98  const Mg3d& mg = *mpaff->get_mg() ;
99  int nz = mg.get_nzone() ; int nzm1 = nz - 1;
100  assert(mg.get_type_r(nzm1) == UNSURR) ;
101  Scalar S_r = *cmp[0] ;
102  if (S_r.get_etat() == ETATZERO) S_r.annule_hard() ;
103  Scalar S_eta = eta() ;
104  if (S_eta.get_etat() == ETATZERO) S_eta.annule_hard() ;
105  S_r.set_spectral_va().ylm() ;
106  S_eta.set_spectral_va().ylm() ;
107  const Base_val& base = S_eta.get_spectral_va().base ;
108  Mtbl_cf sol_part_eta(mg, base) ; sol_part_eta.annule_hard() ;
109  Mtbl_cf sol_part_vr(mg, base) ; sol_part_vr.annule_hard() ;
110  Mtbl_cf solution_hom_un(mg, base) ; solution_hom_un.annule_hard() ;
111  Mtbl_cf solution_hom_deux(mg, base) ; solution_hom_deux.annule_hard() ;
112  Mtbl_cf solution_hom_trois(mg, base) ; solution_hom_trois.annule_hard() ;
113  Mtbl_cf solution_hom_quatre(mg, base) ; solution_hom_quatre.annule_hard() ;
114 
115 
116  // The l_0 component is solved independently // Understand this step
117  //------------------------------------------
118  Scalar sou_l0 = (*cmp[0]) / (1. + lam) ;
119  Param_elliptic param_l0(sou_l0) ;
120  for (int l=0; l<nz; l++)
121  param_l0.set_poisson_vect_r(l, true) ;
122 
123  // Scalar vrl0 = sou_l0.sol_elliptic(param_l0) ;
124  Scalar vrl0 = sou_l0.sol_elliptic_boundary(param_l0, bound_vr, fact_dir, fact_neu) ;
125 
126  // Build-up & inversion of the system for (eta, V^r) in each domain
127  //-----------------------------------------------------------------
128 
129  // Shells
130  //-------
131 
132  int nr ;
133  int nt = mg.get_nt(0) ;
134  int np = mg.get_np(0) ;
135  int l_q = 0 ; int m_q = 0; int base_r = 0 ;
136  double alpha, beta, ech ;
137 
138  assert(num_front+1 < nzm1) ; // Minimum one shell
139  for (int zone=num_front+1 ; zone<nzm1 ; zone++) { //begins the loop on zones
140  nr = mg.get_nr(zone) ;
141  alpha = mpaff->get_alpha()[zone] ;
142  beta = mpaff->get_beta()[zone] ;
143  ech = beta / alpha ;
144  assert (nr > 5) ;
145  assert(nt == mg.get_nt(zone)) ;
146  assert(np == mg.get_np(zone)) ;
147 
148  // Loop on l and m
149  //----------------
150  for (int k=0 ; k<np+1 ; k++) {
151  for (int j=0 ; j<nt ; j++) {
152  base.give_quant_numbers(zone, k, j, m_q, l_q, base_r) ;
153  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
154  int dege1 = 2 ; //degeneracy of eq.1
155  int dege2 = 2 ; //degeneracy of eq.2
156  int nr_eq1 = nr - dege1 ; //Eq.1 is for eta
157  int nr_eq2 = nr - dege2 ; //Eq.2 is for V^r
158  int nrtot = nr_eq1 + nr_eq2 ;
159  Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
160  Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
161  Diff_x2dsdx2 x2d2(base_r, nr); const Matrice& m2d2 = x2d2.get_matrice() ;
162  Diff_xdsdx2 xd2(base_r, nr) ; const Matrice& mxd2 = xd2.get_matrice() ;
163  Diff_dsdx2 d2(base_r, nr) ; const Matrice& md2 = d2.get_matrice() ;
164  Diff_xdsdx xd(base_r, nr) ; const Matrice& mxd = xd.get_matrice() ;
165  Diff_dsdx d1(base_r, nr) ; const Matrice& md = d1.get_matrice() ;
166  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
167 
168  // Building the operator // Which is the eq. from the notes that it is actually implemented?
169  //----------------------
170  for (int lin=0; lin<nr_eq1; lin++) {
171  for (int col=dege1; col<nr; col++)
172  oper.set(lin,col-dege1)
173  = m2d2(lin,col) + 2*ech*mxd2(lin,col) + ech*ech*md2(lin,col)
174  + 2*(mxd(lin,col) + ech*md(lin,col))
175  - (lam+1)*l_q*(l_q+1)*mid(lin,col) ;
176  for (int col=dege2; col<nr; col++)
177  oper.set(lin,col-dege2+nr_eq1)
178  = lam*(mxd(lin,col) + ech*md(lin,col)) + 2*(1+lam)*mid(lin,col) ;
179  }
180  for (int lin=0; lin<nr_eq2; lin++) {
181  for (int col=dege1; col<nr; col++)
182  oper.set(lin+nr_eq1,col-dege1)
183  = -l_q*(l_q+1)*(lam*(mxd(lin,col) + ech*md(lin,col))
184  - (lam+2)*mid(lin,col)) ;
185  for (int col=dege2; col<nr; col++)
186  oper.set(lin+nr_eq1, col-dege2+nr_eq1)
187  = (lam+1)*(m2d2(lin,col) + 2*ech*mxd2(lin,col)
188  + ech*ech*md2(lin,col)
189  + 2*(mxd(lin,col) + ech*md(lin,col)))
190  -(2*(lam+1)+l_q*(l_q+1))*mid(lin,col) ;
191  }
192  oper.set_lu() ;
193 
194  // Filling the r.h.s
195  //------------------
196  Tbl sr(nr) ; sr.set_etat_qcq() ;
197  Tbl seta(nr) ; seta.set_etat_qcq() ;
198  for (int i=0; i<nr; i++) {
199  sr.set(i) = (*S_r.get_spectral_va().c_cf)(zone, k, j, i);
200  seta.set(i) = (*S_eta.get_spectral_va().c_cf)(zone, k, j, i) ;
201  }
202  Tbl xsr= sr ; Tbl x2sr= sr ;
203  Tbl xseta= seta ; Tbl x2seta = seta ;
204  multx2_1d(nr, &x2sr.t, base_r) ; multx_1d(nr, &xsr.t, base_r) ;
205  multx2_1d(nr, &x2seta.t, base_r) ; multx_1d(nr, &xseta.t, base_r) ;
206 
207  for (int i=0; i<nr_eq1; i++)
208  sec_membre.set(i) = alpha*(alpha*x2seta(i) + 2*beta*xseta(i))
209  + beta*beta*seta(i);
210  for (int i=0; i<nr_eq2; i++)
211  sec_membre.set(i+nr_eq1) = beta*beta*sr(i)
212  + alpha*(alpha*x2sr(i) + 2*beta*xsr(i)) ;
213 
214  // Inversion of the "big" operator
215  //--------------------------------
216  Tbl big_res = oper.inverse(sec_membre) ;
217  Tbl res_eta(nr) ; res_eta.set_etat_qcq() ;
218  Tbl res_vr(nr) ; res_vr.set_etat_qcq() ;
219 
220  // Putting coefficients of h and v to individual arrays
221  //-----------------------------------------------------
222  for (int i=0; i<dege1; i++)
223  res_eta.set(i) = 0 ;
224  for (int i=dege1; i<nr; i++)
225  res_eta.set(i) = big_res(i-dege1) ;
226  for (int i=0; i<dege2; i++)
227  res_vr.set(i) = 0 ;
228  for (int i=dege2; i<nr; i++)
229  res_vr.set(i) = big_res(i-dege2+nr_eq1) ;
230 
231  //homogeneous solutions //I do not understand!!!
232  Tbl sol_hom1 = solh(nr, l_q-1, ech, base_r) ;
233  Tbl sol_hom2 = solh(nr, l_q+1, ech, base_r) ;
234  for (int i=0 ; i<nr ; i++) {
235  sol_part_eta.set(zone, k, j, i) = res_eta(i) ;
236  sol_part_vr.set(zone, k, j, i) = res_vr(i) ;
237  solution_hom_un.set(zone, k, j, i) = sol_hom1(0,i) ;
238  solution_hom_deux.set(zone, k, j, i) = sol_hom2(1,i) ;
239  solution_hom_trois.set(zone, k, j, i) = sol_hom2(0,i) ;
240  solution_hom_quatre.set(zone, k, j, i) = sol_hom1(1,i) ;
241  }
242  }
243  }
244  }
245  }
246 
247  // Compactified external domain
248  //-----------------------------
249  nr = mg.get_nr(nzm1) ;
250  assert(nt == mg.get_nt(nzm1)) ;
251  assert(np == mg.get_np(nzm1)) ;
252  alpha = mpaff->get_alpha()[nzm1] ;
253  assert (nr > 4) ;
254  int nr0 = nr - 2 ; //two degrees of freedom less because of division by u^2
255 
256  // Loop on l and m
257  //----------------
258  for (int k=0 ; k<np+1 ; k++) {
259  for (int j=0 ; j<nt ; j++) {
260  base.give_quant_numbers(nzm1, k, j, m_q, l_q, base_r) ;
261  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
262  int dege1 = 1; //degeneracy of eq.1
263  int dege2 = 0; //degeneracy of eq.2
264  int nr_eq1 = nr0 - dege1 ; //Eq.1 is for eta
265  int nr_eq2 = nr0 - dege2 ; //Eq.2 is the div-free condition
266  int nrtot = nr_eq1 + nr_eq2 ;
267  Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
268  Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
269  Diff_x2dsdx2 x2d2(base_r, nr) ; const Matrice& m2d2 = x2d2.get_matrice() ;
270  Diff_xdsdx xd(base_r, nr) ; const Matrice& mxd = xd.get_matrice() ;
271  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
272 
273  // Building the operator
274  //----------------------
275  for (int lin=0; lin<nr_eq1; lin++) {
276  for (int col=dege1; col<nr0; col++)
277  oper.set(lin,col-dege1)
278  = m2d2(lin,col) + 4*mxd(lin,col)
279  + (2-(lam+1)*l_q*(l_q+1))*mid(lin,col) ;
280  for (int col=dege2; col<nr0; col++)
281  oper.set(lin,col-dege2+nr_eq1) =
282  -lam*mxd(lin,col) + 2*mid(lin,col) ;
283  }
284  for (int lin=0; lin<nr_eq2; lin++) {
285  for (int col=dege1; col<nr0; col++)
286  oper.set(lin+nr_eq1,col-dege1)
287  = l_q*(l_q+1)*(lam*mxd(lin,col) + (3*lam+2)*mid(lin,col)) ;
288  for (int col=dege2; col<nr0; col++)
289  oper.set(lin+nr_eq1, col-dege2+nr_eq1)
290  = (lam+1)*(m2d2(lin,col) + 4*mxd(lin,col))
291  - l_q*(l_q+1)*mid(lin,col) ;
292  }
293  oper.set_lu() ;
294 
295  // Filling the r.h.s
296  //------------------
297  for (int i=0; i<nr_eq1; i++)
298  sec_membre.set(i) = (*S_eta.get_spectral_va().c_cf)(nzm1, k, j, i) ;
299  for (int i=0; i<nr_eq2; i++)
300  sec_membre.set(i+nr_eq1) =(*S_r.get_spectral_va().c_cf)(nzm1, k, j, i);
301  Tbl big_res = oper.inverse(sec_membre) ;
302  Tbl res_eta(nr) ; res_eta.set_etat_qcq() ;
303  Tbl res_vr(nr) ; res_vr.set_etat_qcq() ;
304 
305  // Putting coefficients of h and v to individual arrays
306  //-----------------------------------------------------
307  for (int i=0; i<dege1; i++)
308  res_eta.set(i) = 0 ;
309  for (int i=dege1; i<nr0; i++)
310  res_eta.set(i) = big_res(i-dege1) ;
311  res_eta.set(nr0) = 0 ;
312  res_eta.set(nr0+1) = 0 ;
313  for (int i=0; i<dege2; i++)
314  res_vr.set(i) = 0 ;
315  for (int i=dege2; i<nr0; i++)
316  res_vr.set(i) = big_res(i-dege2+nr_eq1) ;
317  res_vr.set(nr0) = 0 ;
318  res_vr.set(nr0+1) = 0 ;
319 
320  // Multiplication by u^2
321  //-----------------------
322  Tbl res_v2(nr) ; res_v2.set_etat_qcq() ;
323  Tbl res_e2(nr) ; res_e2.set_etat_qcq() ;
324  mult2_xm1_1d_cheb(nr, res_eta.t, res_e2.t) ;
325  mult2_xm1_1d_cheb(nr, res_vr.t, res_v2.t) ;
326 
327  // Homogeneous solution (only 1/r^(l+2) and 1/r^l in the CED)
328  Tbl sol_hom1 = solh(nr, l_q-1, 0., base_r) ;
329  Tbl sol_hom2 = solh(nr, l_q+1, 0., base_r) ;
330  for (int i=0 ; i<nr ; i++) {
331  sol_part_eta.set(nzm1, k, j, i) = alpha*alpha*res_e2(i) ;
332  sol_part_vr.set(nzm1, k, j, i) = alpha*alpha*res_v2(i) ;
333  solution_hom_un.set(nzm1, k, j, i) = 0. ;
334  solution_hom_deux.set(nzm1, k, j, i) = sol_hom2(i) ;
335  solution_hom_trois.set(nzm1, k, j, i) = 0. ;
336  solution_hom_quatre.set(nzm1, k, j, i) = sol_hom1(i) ;
337  }
338  }
339  }
340  }
341 
342  // Now let's match everything ...
343  //-------------------------------
344 
345  // Resulting V^r & eta
346  Scalar vr(*mpaff) ; vr.set_etat_qcq() ;
347  vr.set_spectral_base(base) ;
349  Mtbl_cf& cf_vr = *vr.set_spectral_va().c_cf ;
350  cf_vr.annule_hard() ;
351  Scalar het(*mpaff) ; het.set_etat_qcq() ;
352  het.set_spectral_base(base) ;
354  Mtbl_cf& cf_eta = *het.set_spectral_va().c_cf ;
355  cf_eta.annule_hard() ;
356  int taille = 4*(nzm1-num_front)-2 ; //## a verifier
357  Tbl sec_membre(taille) ;
358  Matrice systeme(taille, taille) ;
359  systeme.set_etat_qcq() ;
360  int ligne ; int colonne ;
361 
362  // Loop on l and m
363  //----------------
364  for (int k=0 ; k<np+1 ; k++)
365  for (int j=0 ; j<nt ; j++) {
366  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
367  if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q != 0)) {
368 
369  double f3_eta = lam*l_q + 3*lam + 2 ;
370  double f4_eta = 2 + 2*lam - lam*l_q ;
371  double f3_vr = (l_q+1)*(lam*l_q - 2) ;
372  double f4_vr = l_q*(lam*l_q + lam + 2) ;
373  ligne = 0 ;
374  colonne = 0 ;
375  sec_membre.annule_hard() ;
376  for (int l=0; l<taille; l++)
377  for (int c=0; c<taille; c++)
378  systeme.set(l,c) = 0 ;
379 
380  // First shell
381  nr = mg.get_nr(num_front+1) ;
382  alpha = mpaff->get_alpha()[num_front+1] ;
383  double echelle = mpaff->get_beta()[num_front+1]/alpha ;
384  // Conditions on eta (configuration space)
385  //value and derivative of (x+echelle)^(l-1) at -1
386  systeme.set(ligne, colonne) = pow(echelle-1., double(l_q-1)) ;
387 
388  // value and derivative of 1/(x+echelle) ^(l+2) at -1
389  systeme.set(ligne, colonne+1) = 1/pow(echelle-1., double(l_q+2)) ;
390 
391  //value and derivative of (x+echelle)^(l+1) at -1
392  systeme.set(ligne, colonne+2) = f3_eta*pow(echelle-1., double(l_q+1)) ;
393  // value and derivative of 1/(x+echelle) ^l at -1
394  systeme.set(ligne, colonne+3) = f4_eta/pow(echelle-1., double(l_q)) ;
395  for (int i=0 ; i<nr ; i++)
396  if (i%2 == 0)
397  sec_membre.set(ligne) -= sol_part_eta(num_front+1, k, j, i) ;
398  else sec_membre.set(ligne) += sol_part_eta(num_front+1, k, j, i) ;
399  sec_membre.set(ligne) += bound_eta(num_front+1, k, j, 0) ;
400  ligne++ ;
401 
402  // ... and their couterparts for V^r
403  systeme.set(ligne, colonne) = fact_dir*l_q*pow(echelle-1., double(l_q-1)) + fact_neu*l_q*(l_q-1)*pow(echelle-1., double(l_q-2))/alpha ;
404  systeme.set(ligne, colonne+1) = -fact_dir*(l_q+1)/pow(echelle-1., double(l_q+2)) + fact_neu*(l_q+1)*(l_q+2)/pow(echelle-1., double(l_q+3))/alpha ;
405  systeme.set(ligne, colonne+2) = fact_dir*f3_vr*pow(echelle-1., double(l_q+1)) + fact_neu*f3_vr*(l_q+1)*pow(echelle-1., double(l_q))/alpha ;
406  systeme.set(ligne, colonne+3) = fact_dir*f4_vr/pow(echelle-1., double(l_q)) - fact_neu*(f4_vr*l_q/pow(echelle-1., double(l_q+1)))/alpha ;
407  for (int i=0 ; i<nr ; i++)
408  if (i%2 == 0)
409  sec_membre.set(ligne) -= fact_dir*sol_part_vr(num_front+1, k, j, i) - fact_neu*i*i/alpha*sol_part_vr(num_front+1, k, j, i) ;
410  else sec_membre.set(ligne) += fact_dir*sol_part_vr(num_front+1, k, j, i) - fact_neu*i*i/alpha*sol_part_vr(num_front+1, k, j, i) ;
411  sec_membre.set(ligne) += bound_vr(num_front+1, k, j, 0) ;
412 
413  ligne++ ;
414 
415 
416  // Values at 1
417  // eta
418  //value of (x+echelle)^(l-1) at 1
419  systeme.set(ligne, colonne) = pow(echelle+1., double(l_q-1)) ;
420  // value of 1/(x+echelle) ^(l+2) at 1
421  systeme.set(ligne, colonne+1) = 1./pow(echelle+1., double(l_q+2)) ;
422  //value of (x+echelle)^(l+1) at 1
423  systeme.set(ligne, colonne+2) = f3_eta*pow(echelle+1., double(l_q+1));
424  // value of 1/(x+echelle) ^l at 1
425  systeme.set(ligne, colonne+3) = f4_eta/pow(echelle+1., double(l_q)) ;
426  for (int i=0 ; i<nr ; i++)
427  sec_membre.set(ligne) -= sol_part_eta(num_front+1, k, j, i) ;
428  ligne++ ;
429  // ... and their couterparts for V^r
430  systeme.set(ligne, colonne) = l_q*pow(echelle+1., double(l_q-1)) ;
431  systeme.set(ligne, colonne+1)
432  = -double(l_q+1) / pow(echelle+1., double(l_q+2));
433  systeme.set(ligne, colonne+2) = f3_vr*pow(echelle+1., double(l_q+1)) ;
434  systeme.set(ligne, colonne+3) = f4_vr/pow(echelle+1., double(l_q));
435  for (int i=0 ; i<nr ; i++)
436  sec_membre.set(ligne) -= sol_part_vr(num_front+1, k, j, i) ;
437  ligne++ ;
438 
439  //Derivatives at 1
440  // eta
441  //derivative of (x+echelle)^(l-1) at 1
442  systeme.set(ligne, colonne)
443  = (l_q-1) * pow(echelle+1., double(l_q-2))/alpha ;
444  // derivative of 1/(x+echelle) ^(l+2) at 1
445  systeme.set(ligne, colonne+1)
446  = -(l_q+2) / pow(echelle+1., double(l_q+3))/alpha ;
447  // derivative of (x+echelle)^(l+1) at 1
448  systeme.set(ligne, colonne+2)
449  = f3_eta*(l_q+1) * pow(echelle+1., double(l_q))/alpha;
450  // derivative of 1/(x+echelle) ^l at 1
451  systeme.set(ligne, colonne+3)
452  = -f4_eta*l_q / pow(echelle+1., double(l_q+1))/alpha ;
453  for (int i=0 ; i<nr ; i++)
454  sec_membre.set(ligne) -= i*i/alpha*sol_part_eta(num_front+1, k, j, i) ;
455  ligne++ ;
456  // ... and their couterparts for V^r
457  systeme.set(ligne, colonne)
458  = l_q*(l_q-1) * pow(echelle+1., double(l_q-2))/alpha ;
459  systeme.set(ligne, colonne+1)
460  = (l_q+1)*(l_q+2) / pow(echelle+1., double(l_q+3))/alpha ;
461  systeme.set(ligne, colonne+2)
462  = f3_vr*(l_q+1) * pow(echelle+1., double(l_q))/alpha ;
463  systeme.set(ligne, colonne+3)
464  = -f4_vr*l_q / pow(echelle+1., double(l_q+1))/alpha ;
465  for (int i=0 ; i<nr ; i++)
466  sec_membre.set(ligne) -= i*i/alpha*sol_part_vr(num_front+1, k, j, i) ;
467 
468  colonne += 4 ; // We pass to the next domain
469 
470 
471  // Next shells
472  if (num_front+2<nzm1){
473  for (int zone=num_front+2 ; zone<nzm1 ; zone++) {
474  nr = mg.get_nr(zone) ;
475  alpha = mpaff->get_alpha()[zone] ;
476  echelle = mpaff->get_beta()[zone]/alpha ;
477  ligne -= 3 ;
478  //value of (x+echelle)^(l-1) at -1
479  systeme.set(ligne, colonne) = -pow(echelle-1., double(l_q-1)) ;
480  // value of 1/(x+echelle) ^(l+2) at -1
481  systeme.set(ligne, colonne+1) = -1/pow(echelle-1., double(l_q+2)) ;
482  //value of (x+echelle)^(l+1) at -1
483  systeme.set(ligne, colonne+2) = -f3_eta*pow(echelle-1., double(l_q+1));
484  // value of 1/(x+echelle) ^l at -1
485  systeme.set(ligne, colonne+3) = -f4_eta/pow(echelle-1., double(l_q)) ;
486  for (int i=0 ; i<nr ; i++)
487  if (i%2 == 0)
488  sec_membre.set(ligne) += sol_part_eta(zone, k, j, i) ;
489  else sec_membre.set(ligne) -= sol_part_eta(zone, k, j, i) ;
490  ligne++ ;
491  // ... and their couterparts for V^r
492  systeme.set(ligne, colonne) = -l_q*pow(echelle-1., double(l_q-1)) ;
493  systeme.set(ligne, colonne+1) = (l_q+1)/pow(echelle-1., double(l_q+2));
494  systeme.set(ligne, colonne+2) = -f3_vr*pow(echelle-1., double(l_q+1)) ;
495  systeme.set(ligne, colonne+3) = -f4_vr/pow(echelle-1., double(l_q));
496  for (int i=0 ; i<nr ; i++)
497  if (i%2 == 0)
498  sec_membre.set(ligne) += sol_part_vr(zone, k, j, i) ;
499  else sec_membre.set(ligne) -= sol_part_vr(zone, k, j, i) ;
500  ligne++ ;
501 
502  //derivative of (x+echelle)^(l-1) at -1
503  systeme.set(ligne, colonne)
504  = -(l_q-1)*pow(echelle-1., double(l_q-2))/alpha ;
505  // derivative of 1/(x+echelle) ^(l+2) at -1
506  systeme.set(ligne, colonne+1)
507  = (l_q+2)/pow(echelle-1., double(l_q+3))/alpha ;
508  // derivative of (x+echelle)^(l+1) at -1
509  systeme.set(ligne, colonne+2)
510  = -f3_eta*(l_q+1)*pow(echelle-1., double(l_q))/alpha;
511  // derivative of 1/(x+echelle) ^l at -1
512  systeme.set(ligne, colonne+3)
513  = (f4_eta*l_q/pow(echelle-1., double(l_q+1)))/alpha ;
514  for (int i=0 ; i<nr ; i++)
515  if (i%2 == 0) sec_membre.set(ligne)
516  -= i*i/alpha*sol_part_eta(zone, k, j, i) ;
517  else sec_membre.set(ligne) +=
518  i*i/alpha*sol_part_eta(zone, k, j, i) ;
519  ligne++ ;
520  // ... and their couterparts for V^r
521  systeme.set(ligne, colonne)
522  = -l_q*(l_q-1)*pow(echelle-1., double(l_q-2))/alpha ;
523  systeme.set(ligne, colonne+1)
524  = -(l_q+1)*(l_q+2)/pow(echelle-1., double(l_q+3))/alpha ;
525  systeme.set(ligne, colonne+2)
526  = -f3_vr*(l_q+1)*pow(echelle-1., double(l_q))/alpha ;
527  systeme.set(ligne, colonne+3)
528  = (f4_vr*l_q/pow(echelle-1., double(l_q+1)))/alpha ;
529  for (int i=0 ; i<nr ; i++)
530  if (i%2 == 0) sec_membre.set(ligne)
531  -= i*i/alpha*sol_part_vr(zone, k, j, i) ;
532  else sec_membre.set(ligne) +=
533  i*i/alpha*sol_part_vr(zone, k, j, i) ;
534  ligne++ ;
535 
536  //value of (x+echelle)^(l-1) at 1
537  systeme.set(ligne, colonne) = pow(echelle+1., double(l_q-1)) ;
538  // value of 1/(x+echelle) ^(l+2) at 1
539  systeme.set(ligne, colonne+1) = 1./pow(echelle+1., double(l_q+2)) ;
540  //value of (x+echelle)^(l+1) at 1
541  systeme.set(ligne, colonne+2) = f3_eta*pow(echelle+1., double(l_q+1));
542  // value of 1/(x+echelle) ^l at 1
543  systeme.set(ligne, colonne+3) = f4_eta/pow(echelle+1., double(l_q)) ;
544  for (int i=0 ; i<nr ; i++)
545  sec_membre.set(ligne) -= sol_part_eta(zone, k, j, i) ;
546  ligne++ ;
547  // ... and their couterparts for V^r
548  systeme.set(ligne, colonne) = l_q*pow(echelle+1., double(l_q-1)) ;
549  systeme.set(ligne, colonne+1)
550  = -double(l_q+1) / pow(echelle+1., double(l_q+2));
551  systeme.set(ligne, colonne+2) = f3_vr*pow(echelle+1., double(l_q+1)) ;
552  systeme.set(ligne, colonne+3) = f4_vr/pow(echelle+1., double(l_q));
553  for (int i=0 ; i<nr ; i++)
554  sec_membre.set(ligne) -= sol_part_vr(zone, k, j, i) ;
555  ligne++ ;
556 
557  //derivative of (x+echelle)^(l-1) at 1
558  systeme.set(ligne, colonne)
559  = (l_q-1) * pow(echelle+1., double(l_q-2))/alpha ;
560  // derivative of 1/(x+echelle) ^(l+2) at 1
561  systeme.set(ligne, colonne+1)
562  = -(l_q+2) / pow(echelle+1., double(l_q+3))/alpha ;
563  // derivative of (x+echelle)^(l+1) at 1
564  systeme.set(ligne, colonne+2)
565  = f3_eta*(l_q+1) * pow(echelle+1., double(l_q))/alpha;
566  // derivative of 1/(x+echelle) ^l at 1
567  systeme.set(ligne, colonne+3)
568  = -f4_eta*l_q / pow(echelle+1., double(l_q+1))/alpha ;
569  for (int i=0 ; i<nr ; i++)
570  sec_membre.set(ligne) -= i*i/alpha*sol_part_eta(zone, k, j, i) ;
571  ligne++ ;
572  // ... and their couterparts for V^r
573  systeme.set(ligne, colonne)
574  = l_q*(l_q-1) * pow(echelle+1., double(l_q-2))/alpha ;
575  systeme.set(ligne, colonne+1)
576  = (l_q+1)*(l_q+2) / pow(echelle+1., double(l_q+3))/alpha ;
577  systeme.set(ligne, colonne+2)
578  = f3_vr*(l_q+1) * pow(echelle+1., double(l_q))/alpha ;
579  systeme.set(ligne, colonne+3)
580  = -f4_vr*l_q / pow(echelle+1., double(l_q+1))/alpha ;
581  for (int i=0 ; i<nr ; i++)
582  sec_membre.set(ligne) -= i*i/alpha*sol_part_vr(zone, k, j, i) ;
583 
584  colonne += 4 ;
585  }
586  }
587  //Compactified external domain
588  nr = mg.get_nr(nzm1) ;
589 
590  alpha = mpaff->get_alpha()[nzm1] ;
591  ligne -= 3 ;
592  //value of (x-1)^(l+2) at -1 :
593  systeme.set(ligne, colonne) = -pow(-2, double(l_q+2)) ;
594  //value of (x-1)^l at -1 :
595  systeme.set(ligne, colonne+1) = -f4_eta*pow(-2, double(l_q)) ;
596  for (int i=0 ; i<nr ; i++)
597  if (i%2 == 0) sec_membre.set(ligne) += sol_part_eta(nzm1, k, j, i) ;
598  else sec_membre.set(ligne) -= sol_part_eta(nzm1, k, j, i) ;
599  //... and of its couterpart for V^r
600  systeme.set(ligne+1, colonne) = double(l_q+1)*pow(-2, double(l_q+2)) ;
601  systeme.set(ligne+1, colonne+1) = -f4_vr*pow(-2, double(l_q)) ;
602  for (int i=0 ; i<nr ; i++)
603  if (i%2 == 0) sec_membre.set(ligne+1) += sol_part_vr(nzm1, k, j, i) ;
604  else sec_membre.set(ligne+1) -= sol_part_vr(nzm1, k, j, i) ;
605 
606  ligne += 2 ;
607  //derivative of (x-1)^(l+2) at -1 :
608  systeme.set(ligne, colonne) = alpha*(l_q+2)*pow(-2, double(l_q+3)) ;
609  //derivative of (x-1)^l at -1 :
610  systeme.set(ligne, colonne+1) = alpha*l_q*f4_eta*pow(-2, double(l_q+1)) ;
611  for (int i=0 ; i<nr ; i++)
612  if (i%2 == 0) sec_membre.set(ligne)
613  -= -4*alpha*i*i*sol_part_eta(nzm1, k, j, i) ;
614  else sec_membre.set(ligne)
615  += -4*alpha*i*i*sol_part_eta(nzm1, k, j, i) ;
616  //... and of its couterpart for V^r
617  systeme.set(ligne+1, colonne)
618  = -alpha*double((l_q+1)*(l_q+2))*pow(-2, double(l_q+3)) ;
619  systeme.set(ligne+1, colonne+1)
620  = alpha*double(l_q)*f4_vr*pow(-2, double(l_q+1)) ;
621  for (int i=0 ; i<nr ; i++)
622  if (i%2 == 0) sec_membre.set(ligne+1)
623  -= -4*alpha*i*i*sol_part_vr(nzm1, k, j, i) ;
624  else sec_membre.set(ligne+1)
625  += -4*alpha*i*i*sol_part_vr(nzm1, k, j, i) ;
626 
627  // Solution of the system giving the coefficients for the homogeneous
628  // solutions
629  //-------------------------------------------------------------------
630  if (taille > 2) systeme.set_band(5,5) ;
631  else systeme.set_band(1,1) ;
632  systeme.set_lu() ;
633  Tbl facteurs(systeme.inverse(sec_membre)) ;
634  int conte = 0 ;
635 
636  // everything is put to the right place, the same combination of hom.
637  // solutions (with some l or -(l+1) factors) must be used for V^r
638  //-------------------------------------------------------------------
639 
640  for (int zone=1 ; zone<nzm1 ; zone++) { //shells
641  nr = mg.get_nr(zone) ;
642  for (int i=0 ; i<nr ; i++) {
643  cf_eta.set(zone, k, j, i) =
644  sol_part_eta(zone, k, j, i)
645  +facteurs(conte)*solution_hom_un(zone, k, j, i)
646  +facteurs(conte+1)*solution_hom_deux(zone, k, j, i)
647  +facteurs(conte+2)*f3_eta*solution_hom_trois(zone, k, j, i)
648  +facteurs(conte+3)*f4_eta*solution_hom_quatre(zone, k, j, i) ;
649  cf_vr.set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
650  +double(l_q)*facteurs(conte)*solution_hom_un(zone, k, j, i)
651  -double(l_q+1)*facteurs(conte+1)*solution_hom_deux(zone, k, j, i) // What is the origin of these factors?!
652  +f3_vr*facteurs(conte+2)*solution_hom_trois(zone, k, j, i)
653  +f4_vr*facteurs(conte+3)*solution_hom_quatre(zone, k, j, i) ;
654  }
655  conte+=4 ;
656  }
657  nr = mg.get_nr(nz-1) ; //compactified external domain
658  for (int i=0 ; i<nr ; i++) {
659  cf_eta.set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
660  +facteurs(conte)*solution_hom_deux(nzm1, k, j, i)
661  +f4_eta*facteurs(conte+1)*solution_hom_quatre(nzm1, k, j, i) ;
662  cf_vr.set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
663  -double(l_q+1)*facteurs(conte)*solution_hom_deux(nzm1, k, j, i)
664  +f4_vr*facteurs(conte+1)*solution_hom_quatre(nzm1, k, j, i) ;
665 
666  }
667  } // End of nullite_plm
668  } //End of loop on theta
669  vr.set_spectral_va().ylm_i() ;
670  vr += vrl0 ;
671  het.set_spectral_va().ylm_i() ;
672 
673  Valeur temp_mu(mg.get_angu()) ;
674  temp_mu = bound_mu ;
675  const Valeur& limit_mu (temp_mu) ;
676 
677  resu.set_vr_eta_mu(vr, 0*het, mu().poisson_dirichlet(limit_mu, num_front)) ;
678 
679  return ;
680 }
681 
682 
683 Vector Vector::poisson_dirichlet(double lam, const Valeur& bound_vr,
684  const Valeur& bound_vt, const Valeur& bound_vp,
685  int num_front) const {
686 
687  Vector resu(*mp, CON, triad) ;
688  resu = poisson_robin(lam, bound_vr, bound_vt, bound_vp, 1., 0., num_front) ;
689 
690  return resu ;
691 
692 }
693 
694 Vector Vector::poisson_neumann(double lam, const Valeur& bound_vr,
695  const Valeur& bound_vt, const Valeur& bound_vp,
696  int num_front) const {
697 
698  Vector resu(*mp, CON, triad) ;
699  resu = poisson_robin(lam, bound_vr, bound_vt, bound_vp, 0., 1., num_front) ;
700 
701  return resu ;
702 
703 }
704 
705 Vector Vector::poisson_robin(double lam, const Valeur& bound_vr,
706  const Valeur& bound_vt, const Valeur& bound_vp,
707  double fact_dir, double fact_neu,
708  int num_front) const {
709 
710 
711  // Boundary condition for V^r //Construction of a Mtbl_cf from Valeur with Ylm coefficients
712  Valeur limit_vr (bound_vr) ;
713 
714  limit_vr.coef() ;
715  limit_vr.ylm() ; // Spherical harmonics transform.
716  Mtbl_cf lim_vr (*(limit_vr.c_cf)) ;
717 
718  // bound_vt and bound_vp are only known at the boundary --> we fill
719  // all the zones extending the values at the boundary before calling to poisson_angu.
720  Scalar temp_vt (*mp) ;
721  Scalar temp_vp (*mp) ;
722  temp_vt.annule_hard() ;
723  temp_vp.annule_hard() ;
724  int nz = mp->get_mg()->get_nzone() ;
725  for (int l=0; l<nz; l++)
726  for (int j=0; j<mp->get_mg()->get_nt(l); j++)
727  for (int k=0; k<mp->get_mg()->get_np(l); k++) {
728  temp_vt.set_grid_point(l, k, j, 0) = bound_vt(1, k, j, 0) ;
729  temp_vp.set_grid_point(l, k, j, 0) = bound_vp(1, k, j, 0) ;
730  }
731  temp_vt.set_spectral_va().set_base(bound_vt.base) ; // We set the basis
732  temp_vp.set_spectral_va().set_base(bound_vp.base) ;
733 
734  cout << "temp_vp" << endl << norme(temp_vp) << endl ;
735 
736  //Source for eta
737  Scalar source_eta(*mp) ;
738  Scalar vtstant (temp_vt) ;
739  vtstant.div_tant() ;
740  source_eta = temp_vt.dsdt() + vtstant + temp_vp.stdsdp() ;
741 
742  //Source for mu
743  Scalar source_mu(*mp) ;
744  Scalar vpstant (temp_vp) ;
745  vpstant.div_tant() ;
746  source_mu = temp_vp.dsdt() + vpstant - temp_vt.stdsdp() ; //There was a wrong sign here
747 
748  Scalar temp_mu (source_mu.poisson_angu()) ;
749  Scalar temp_eta (source_eta.poisson_angu()) ;
750 
751  // Boundary condition for mu
752  Valeur limit_mu ((*mp).get_mg()->get_angu() ) ;
753  int nnp = (*mp).get_mg()->get_np(1) ;
754  int nnt = (*mp).get_mg()->get_nt(1) ;
755  limit_mu= 1 ;
756  for (int k=0 ; k<nnp ; k++)
757  for (int j=0 ; j<nnt ; j++)
758  limit_mu.set(1, k, j, 0) = temp_mu.val_grid_point(1, k, j, 0) ;
759  limit_mu.set_base(temp_mu.get_spectral_va().get_base()) ;
760 
761  limit_mu.coef() ;
762  limit_mu.ylm() ; // Spherical harmonics transform.
763  Mtbl_cf lim_mu (*(limit_mu.c_cf)) ;
764 
765  // Boundary condition for eta
766  Valeur limit_eta ((*mp).get_mg()->get_angu() ) ;
767  limit_eta = 1 ;
768  for (int k=0 ; k<nnp ; k++)
769  for (int j=0 ; j<nnt ; j++)
770  limit_eta.set(1, k, j, 0) = temp_eta.val_grid_point(1, k, j, 0) ;
771  limit_eta.set_base(temp_eta.get_spectral_va().get_base()) ;
772 
773  limit_eta.coef() ;
774  limit_eta.ylm() ; // Spherical harmonics transform.
775  Mtbl_cf lim_eta (*(limit_eta.c_cf)) ;
776 
777 
778  // Call to poisson_boundary(...)
779  bool nullite = true ;
780  for (int i=0; i<3; i++) {
781  assert(cmp[i]->check_dzpuis(4)) ;
782  if (cmp[i]->get_etat() != ETATZERO || bound_vr.get_etat() != ETATZERO ||
783  bound_vt.get_etat() != ETATZERO || bound_vp.get_etat() != ETATZERO)
784  nullite = false ;
785  }
786 
787  Vector resu(*mp, CON, triad) ;
788  if (nullite)
789  resu.set_etat_zero() ;
790  else
791  poisson_boundary(lam, lim_vr, lim_eta, lim_mu, num_front, fact_dir,
792  fact_neu, resu) ;
793 
794  return resu ;
795 
796 }
797 }
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_x2dsdx2.C:98
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
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:607
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:134
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
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
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:151
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:395
Lorene prototypes.
Definition: app_hor.h:67
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:304
void poisson_boundary(double lambda, const Mtbl_cf &limit_vr, const Mtbl_cf &limit_eta, const Mtbl_cf &limit_mu, int num_front, double fact_dir, double fact_neu, Vector &resu) const
Solves the vector Poisson equation with *this as a source with a boundary condition on the excised sp...
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition: matrice.C:427
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:141
const Scalar & dsdt() const
Returns of *this .
Definition: scalar_deriv.C:208
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_dsdx.C:97
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Vector poisson_robin(double lambda, const Valeur &limit_vr, const Valeur &limit_vt, const Valeur &limit_vp, double fact_dir, double fact_neu, int num_front) const
Solves the vector Poisson equation with *this as a source with a boundary condition on the excised sp...
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tensor.h:309
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:129
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:359
Tensor field of valence 1.
Definition: vector.h:188
Scalar poisson_angu(double lambda=0) const
Solves the (generalized) angular Poisson equation with *this as source.
Definition: scalar_pde.C:203
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition: valeur.C:813
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:386
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
Definition: valeur.h:763
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
Class for the elementary differential operator Identity (see the base class Diff ).
Definition: diff.h:210
int get_etat() const
Returns the logical state.
Definition: valeur.h:760
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_xdsdx2.C:100
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:490
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:484
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written: .
Definition: vector_etamu.C:72
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition: mg3d.C:604
Matrix handling.
Definition: matrice.h:152
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:315
void set_vr_eta_mu(const Scalar &vr_i, const Scalar &eta_i, const Scalar &mu_i)
Defines the components through potentials and (see members p_eta and p_mu ), as well as the compon...
Definition: vector_etamu.C:198
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition: tensor.h:321
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_dsdx2.C:94
double * t
The array of double.
Definition: tbl.h:176
const double * get_beta() const
Returns the pointer on the array beta.
Definition: map_af.C:611
Scalar sol_elliptic_boundary(Param_elliptic &params, const Mtbl_cf &bound, double fact_dir, double fact_neu) const
Resolution of a general elliptic equation, putting zero at infinity and with inner boundary condition...
Definition: scalar_pde.C:259
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
const Scalar & stdsdp() const
Returns of *this .
Definition: scalar_deriv.C:238
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
This class contains the parameters needed to call the general elliptic solver.
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:367
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:803
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Multi-domain grid.
Definition: grilles.h:279
Bases of the spectral expansions.
Definition: base_val.h:325
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition: scalar.h:690
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition: mtbl_cf.C:315
Affine radial mapping.
Definition: map.h:2048
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:409
Vector poisson_neumann(double lambda, const Valeur &limit_vr, const Valeur &limit_vt, const Valeur &limit_vp, int num_front) const
Solves the vector Poisson equation with *this as a source with a boundary condition on the excised sp...
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:196
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written: .
Definition: vector_etamu.C:107
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:178
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:531
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:506
Basic array class.
Definition: tbl.h:164
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:491
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:301
void div_tant()
Division by .
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:375
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_xdsdx.C:101
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607
void set_poisson_vect_r(int zone, bool only_l_zero=false)
Sets the operator to in all domains, for ; and to in all domains otherwise.
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition: valeur.h:373