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