LORENE
vector_poisson_boundary2.C
1 /*
2  * Method for vector Poisson equation inverting eqs. for V^r and eta as a block.
3  *
4  * (see file vector.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2005 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_poisson_boundary2_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/vector_poisson_boundary2.C,v 1.9 2014/10/13 08:53:45 j_novak Exp $" ;
29 
30 /*
31  * $Id: vector_poisson_boundary2.C,v 1.9 2014/10/13 08:53:45 j_novak Exp $
32  * $Log: vector_poisson_boundary2.C,v $
33  * Revision 1.9 2014/10/13 08:53:45 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.8 2014/10/06 15:13:21 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.7 2008/08/20 15:07:36 n_vasset
40  * Cleaning up the code...
41  *
42  * Revision 1.5 2007/09/05 12:35:18 j_novak
43  * Homogeneous solutions are no longer obtained through the analytic formula, but
44  * by solving (again) the operator system with almost zero as r.h.s. This is to
45  * lower the condition number of the matching system.
46  *
47  * Revision 1.4 2007/01/23 17:08:46 j_novak
48  * New function pois_vect_r0.C to solve the l=0 part of the vector Poisson
49  * equation, which involves only the r-component.
50  *
51  * Revision 1.3 2005/12/30 13:39:38 j_novak
52  * Changed the Galerkin base in the nucleus to (hopefully) stabilise the solver
53  * when used in an iteration. Similar changes in the CED too.
54  *
55  * Revision 1.2 2005/02/15 15:43:18 j_novak
56  * First version of the block inversion for the vector Poisson equation (method 6).
57  *
58  *
59  * $Header: /cvsroot/Lorene/C++/Source/Tensor/vector_poisson_boundary2.C,v 1.9 2014/10/13 08:53:45 j_novak Exp $
60  *
61  */
62 
63 // C headers
64 #include <cassert>
65 #include <cstdlib>
66 #include <cmath>
67 
68 // Lorene headers
69 #include "metric.h"
70 #include "diff.h"
71 #include "param_elliptic.h"
72 #include "proto.h"
73 #include "graphique.h"
74 #include "map.h"
75 #include "utilitaires.h"
76 
77 
78 
79 namespace Lorene {
80 void Vector::poisson_boundary2(double lam, Vector& resu, Scalar boundvr, Scalar boundeta, Scalar boundmu, double dir_vr, double neum_vr, double dir_eta, double neum_eta, double dir_mu, double neum_mu) const {
81 
82  const Map_af* mpaff = dynamic_cast<const Map_af*>(mp) ;
83 #ifndef NDEBUG
84  for (int i=0; i<3; i++)
85  assert(cmp[i]->check_dzpuis(4)) ;
86  // All this has a meaning only for spherical components:
87  const Base_vect_spher* bvs = dynamic_cast<const Base_vect_spher*>(triad) ;
88  assert(bvs != 0x0) ;
89  //## ... and affine mapping, for the moment!
90  assert( mpaff != 0x0) ;
91 #endif
92 
93  if (fabs(lam + 1.) < 1.e-8) {
94 
95  cout <<" lambda = -1 is not supported by this code!" << endl;
96  abort();
97  }
98 
99 
100 
101  // Extraction of Mtbl_cf objects from boundary informations;
102  Scalar bound_vr2 = boundvr;
103  bound_vr2.set_spectral_va().ylm();
104  Scalar bound_eta2 = boundeta;
105  bound_eta2.set_spectral_va().ylm();
106  Scalar bound_mu2 = boundmu; bound_mu2.set_spectral_va().ylm();
107 
108  const Mtbl_cf *bound_vr = bound_vr2.get_spectral_va().c_cf;
109  const Mtbl_cf *bound_mu = bound_mu2.get_spectral_va().c_cf;
110  const Mtbl_cf *bound_eta = bound_eta2.get_spectral_va().c_cf;
111 
112 
113 
114 
115  // Some working objects
116  //---------------------
117  const Mg3d& mg = *mpaff->get_mg() ;
118  int nz = mg.get_nzone() ; int nzm1 = nz - 1;
119  assert(mg.get_type_r(0) == RARE) ;
120  assert(mg.get_type_r(nzm1) == UNSURR) ;
121  Scalar S_r = *cmp[0] ;
122  Scalar S_eta = eta() ;
123  Scalar het(*mpaff) ;
124  Scalar vr(*mpaff) ;
125 
126  if (S_r.get_etat() == ETATZERO) {
127  if (S_eta.get_etat() == ETATZERO) {
128 
129  S_r.annule_hard() ;
130  S_r.set_spectral_base(S_eta.get_spectral_base()) ;
131  S_eta.annule_hard() ;
132  S_eta.set_spectral_base(S_r.get_spectral_base()) ;
133  }
134  else {
135  S_r.annule_hard() ;
136  S_r.set_spectral_base(S_eta.get_spectral_base()) ;
137  }
138  }
139  if (S_eta.get_etat() == ETATZERO) {
140  S_eta.annule_hard() ;
141  S_eta.set_spectral_base(S_r.get_spectral_base()) ;
142  }
143 
144  S_r.set_spectral_va().ylm() ;
145  S_eta.set_spectral_va().ylm() ;
146  const Base_val& base = S_eta.get_spectral_va().base ;
147  Mtbl_cf sol_part_eta(mg, base) ; sol_part_eta.annule_hard() ;
148  Mtbl_cf sol_part_vr(mg, base) ; sol_part_vr.annule_hard() ;
149  Mtbl_cf sol_hom_un_eta(mg, base) ; sol_hom_un_eta.annule_hard() ;
150  Mtbl_cf sol_hom_deux_eta(mg, base) ; sol_hom_deux_eta.annule_hard() ;
151  Mtbl_cf sol_hom_trois_eta(mg, base) ; sol_hom_trois_eta.annule_hard() ;
152  Mtbl_cf sol_hom_quatre_eta(mg, base) ; sol_hom_quatre_eta.annule_hard() ;
153  Mtbl_cf sol_hom_un_vr(mg, base) ; sol_hom_un_vr.annule_hard() ;
154  Mtbl_cf sol_hom_deux_vr(mg, base) ; sol_hom_deux_vr.annule_hard() ;
155  Mtbl_cf sol_hom_trois_vr(mg, base) ; sol_hom_trois_vr.annule_hard() ;
156  Mtbl_cf sol_hom_quatre_vr(mg, base) ; sol_hom_quatre_vr.annule_hard() ;
157 
158  // Solution of the l=0 part for Vr
159  //---------------------------------
160  Scalar sou_l0 = (*cmp[0]) / (1. + lam) ;
161  sou_l0.set_spectral_va().ylm() ;
162 
163  Param_elliptic param_l0(sou_l0) ;
164  for (int l=0; l<nz; l++){
165  param_l0.set_poisson_vect_r(l, true) ;
166  }
167 
168 
169  Scalar vrl0 = sou_l0.sol_elliptic_boundary(param_l0, *bound_vr, dir_vr, neum_vr) ;
170 
171 
172  // Build-up & inversion of the system for (eta, V^r) in each domain (except for the nucleus!)
173  //-----------------------------------------------------------------
174 
175 
176  // Shells
177  //-------
178 
179  int nr = mg.get_nr(1) ;
180  int nt = mg.get_nt(1) ;
181  int np = mg.get_np(1) ;
182  double alpha = mpaff->get_alpha()[1] ; double alp2 = alpha*alpha ;
183  double beta = mpaff->get_beta()[1] ;
184 
185  int l_q = 0 ; int m_q = 0; int base_r = 0 ;
186  for (int zone=1 ; zone<nzm1 ; zone++) {
187  nr = mg.get_nr(zone) ;
188  assert (nr > 5) ;
189  assert(nt == mg.get_nt(zone)) ;
190  assert(np == mg.get_np(zone)) ;
191  alpha = mpaff->get_alpha()[zone] ;
192  beta = mpaff->get_beta()[zone] ;
193  double ech = beta / alpha ;
194 
195  // Loop on l and m
196  //----------------
197  for (int k=0 ; k<np+1 ; k++) {
198  for (int j=0 ; j<nt ; j++) {
199  base.give_quant_numbers(zone, k, j, m_q, l_q, base_r) ;
200  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
201  int dege1 = 2 ; //degeneracy of eq.1
202  int dege2 = 2 ; //degeneracy of eq.2
203  int nr_eq1 = nr - dege1 ; //Eq.1 is for eta
204  int nr_eq2 = nr - dege2 ; //Eq.2 is for V^r
205  int nrtot = 2*nr ;
206  Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
207  Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
208  Diff_x2dsdx2 x2d2(base_r, nr); const Matrice& m2d2 = x2d2.get_matrice() ;
209  Diff_xdsdx2 xd2(base_r, nr) ; const Matrice& mxd2 = xd2.get_matrice() ;
210  Diff_dsdx2 d2(base_r, nr) ; const Matrice& md2 = d2.get_matrice() ;
211  Diff_xdsdx xd(base_r, nr) ; const Matrice& mxd = xd.get_matrice() ;
212  Diff_dsdx d1(base_r, nr) ; const Matrice& md = d1.get_matrice() ;
213  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
214 
215  // Building the operator
216  //----------------------
217  for (int lin=0; lin<nr_eq1; lin++) {
218  for (int col=0; col<nr; col++)
219  oper.set(lin,col)
220  = m2d2(lin,col) + 2*ech*mxd2(lin,col) + ech*ech*md2(lin,col)
221  + 2*(mxd(lin,col) + ech*md(lin,col))
222  - (lam+1)*l_q*(l_q+1)*mid(lin,col) ;
223  for (int col=0; col<nr; col++)
224  oper.set(lin,col+nr)
225  = lam*(mxd(lin,col) + ech*md(lin,col))
226  + 2*(1+lam)*mid(lin,col) ;
227  }
228  oper.set(nr_eq1, 0) = 1 ;
229  for (int col=1; col<2*nr; col++)
230  oper.set(nr_eq1, col) = 0 ;
231  oper.set(nr_eq1+1, 0) = 0 ;
232  oper.set(nr_eq1+1, 1) = 1 ;
233  for (int col=2; col<2*nr; col++)
234  oper.set(nr_eq1+1, col) = 0 ;
235  for (int lin=0; lin<nr_eq2; lin++) {
236  for (int col=0; col<nr; col++)
237  oper.set(lin+nr,col)
238  = -l_q*(l_q+1)*(lam*(mxd(lin,col) + ech*md(lin,col))
239  - (lam+2)*mid(lin,col)) ;
240  for (int col=0; col<nr; col++)
241  oper.set(lin+nr, col+nr)
242  = (lam+1)*(m2d2(lin,col) + 2*ech*mxd2(lin,col)
243  + ech*ech*md2(lin,col)
244  + 2*(mxd(lin,col) + ech*md(lin,col)))
245  -(2*(lam+1)+l_q*(l_q+1))*mid(lin,col) ;
246  }
247  for (int col=0; col<nr; col++)
248  oper.set(nr+nr_eq2, col) = 0 ;
249  oper.set(nr+nr_eq2, nr) = 1 ;
250  for (int col=nr+1; col<2*nr; col++)
251  oper.set(nr+nr_eq2, col) = 0 ;
252  for (int col=0; col<nr+1; col++)
253  oper.set(nr+nr_eq2+1, col) = 0 ;
254  oper.set(nr+nr_eq2+1, nr+1) = 1 ;
255  for (int col=nr+2; col<2*nr; col++)
256  oper.set(nr+nr_eq2+1, col) = 0 ;
257 
258  oper.set_lu() ;
259 
260  // Filling the r.h.s
261  //------------------
262  Tbl sr(nr) ; sr.set_etat_qcq() ;
263  Tbl seta(nr) ; seta.set_etat_qcq() ;
264  for (int i=0; i<nr; i++) {
265  sr.set(i) = (*S_r.get_spectral_va().c_cf)(zone, k, j, i);
266  seta.set(i) = (*S_eta.get_spectral_va().c_cf)(zone, k, j, i) ;
267  }
268  Tbl xsr= sr ; Tbl x2sr= sr ;
269  Tbl xseta= seta ; Tbl x2seta = seta ;
270  multx2_1d(nr, &x2sr.t, base_r) ; multx_1d(nr, &xsr.t, base_r) ;
271  multx2_1d(nr, &x2seta.t, base_r) ; multx_1d(nr, &xseta.t, base_r) ;
272 
273  for (int i=0; i<nr_eq1; i++)
274  sec_membre.set(i) = alpha*(alpha*x2seta(i) + 2*beta*xseta(i))
275  + beta*beta*seta(i);
276  sec_membre.set(nr_eq1) = 0 ;
277  sec_membre.set(nr_eq1+1) = 0 ;
278  for (int i=0; i<nr_eq2; i++)
279  sec_membre.set(i+nr) = beta*beta*sr(i)
280  + alpha*(alpha*x2sr(i) + 2*beta*xsr(i)) ;
281  sec_membre.set(nr+nr_eq2) = 0 ;
282  sec_membre.set(nr+nr_eq2+1) = 0 ;
283 
284  // Inversion of the "big" operator
285  //--------------------------------
286  Tbl big_res = oper.inverse(sec_membre) ;
287  for (int i=0; i<nr; i++) {
288  sol_part_eta.set(zone, k, j, i) = big_res(i) ;
289  sol_part_vr.set(zone, k, j, i) = big_res(nr+i) ;
290  }
291 
292  // Getting the homogeneous solutions
293  //----------------------------------
294  sec_membre.annule_hard() ;
295  sec_membre.set(nr_eq1) = 1 ;
296  big_res = oper.inverse(sec_membre) ;
297  for (int i=0 ; i<nr ; i++) {
298  sol_hom_un_eta.set(zone, k, j, i) = big_res(i) ;
299  sol_hom_un_vr.set(zone, k, j, i) = big_res(nr+i) ;
300  }
301  sec_membre.set(nr_eq1) = 0 ;
302  sec_membre.set(nr_eq1+1) = 1 ;
303  big_res = oper.inverse(sec_membre) ;
304  for (int i=0 ; i<nr ; i++) {
305  sol_hom_deux_eta.set(zone, k, j, i) = big_res(i) ;
306  sol_hom_deux_vr.set(zone, k, j, i) = big_res(nr+i) ;
307  }
308  sec_membre.set(nr_eq1+1) = 0 ;
309  sec_membre.set(nr+nr_eq2) = 1 ;
310  big_res = oper.inverse(sec_membre) ;
311  for (int i=0 ; i<nr ; i++) {
312  sol_hom_trois_eta.set(zone, k, j, i) = big_res(i) ;
313  sol_hom_trois_vr.set(zone, k, j, i) = big_res(nr+i) ;
314  }
315  sec_membre.set(nr+nr_eq2) = 0 ;
316  sec_membre.set(nr+nr_eq2+1) = 1 ;
317  big_res = oper.inverse(sec_membre) ;
318  for (int i=0 ; i<nr ; i++) {
319  sol_hom_quatre_eta.set(zone, k, j, i) = big_res(i) ;
320  sol_hom_quatre_vr.set(zone, k, j, i) = big_res(nr+i) ;
321  }
322 
323  }
324  }
325  }
326  }
327 
328  // Compactified external domain
329  //-----------------------------
330  nr = mg.get_nr(nzm1) ;
331  assert(nt == mg.get_nt(nzm1)) ;
332  assert(np == mg.get_np(nzm1)) ;
333  alpha = mpaff->get_alpha()[nzm1] ; alp2 = alpha*alpha ;
334  assert (nr > 4) ;
335 
336  // Loop on l and m
337  //----------------
338  for (int k=0 ; k<np+1 ; k++) {
339  for (int j=0 ; j<nt ; j++) {
340  base.give_quant_numbers(nzm1, k, j, m_q, l_q, base_r) ;
341  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
342  int dege1 = 3; //degeneracy of eq.1
343  int dege2 = (l_q == 1 ? 2 : 3); //degeneracy of eq.2
344  int nr_eq1 = nr - dege1 ; //Eq.1 is for eta
345  int nr_eq2 = nr - dege2 ; //Eq.2 is the div-free condition
346  int nrtot = 2*nr ;
347  Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
348  Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
349  Diff_dsdx2 d2(base_r, nr) ; const Matrice& md2 = d2.get_matrice() ;
350  Diff_sxdsdx sxd(base_r, nr) ; const Matrice& mxd = sxd.get_matrice() ;
351  Diff_sx2 sx2(base_r, nr) ; const Matrice& ms2 = sx2.get_matrice() ;
352 
353  // Building the operator
354  //----------------------
355  for (int lin=0; lin<nr_eq1; lin++) {
356  for (int col=0; col<nr; col++)
357  oper.set(lin,col)
358  = (md2(lin,col) - (lam+1)*l_q*(l_q+1)*ms2(lin,col))/alp2 ;
359  for (int col=0; col<nr; col++)
360  oper.set(lin,col+nr) =
361  (-lam*mxd(lin,col) + 2*(1+lam)*ms2(lin,col)) / alp2 ;
362  }
363  for (int col=0; col<nr; col++)
364  oper.set(nr_eq1, col) = 1 ;
365  for (int col=nr; col<nrtot; col++)
366  oper.set(nr_eq1, col) = 0 ;
367  int pari = -1 ;
368  for (int col=0; col<nr; col++) {
369  oper.set(nr_eq1+1, col) = pari*col*col ;
370  pari = pari ;
371  }
372  for (int col=nr; col<nrtot; col++)
373  oper.set(nr_eq1+1, col) = 0 ;
374  oper.set(nr_eq1+2, 0) = 1 ;
375  for (int col=1; col<nrtot; col++)
376  oper.set(nr_eq1+2, col) = 0 ;
377  for (int lin=0; lin<nr_eq2; lin++) {
378  for (int col=0; col<nr; col++)
379  oper.set(lin+nr,col) = (l_q*(l_q+1)*(lam*mxd(lin,col)
380  + (lam+2)*ms2(lin,col))) / alp2 ;
381  for (int col=0; col<nr; col++)
382  oper.set(lin+nr, col+nr) = ((lam+1)*md2(lin,col)
383  - (2*(lam+1) + l_q*(l_q+1))*ms2(lin,col)) / alp2 ;
384  }
385  for (int col=0; col<nr; col++)
386  oper.set(nr+nr_eq2, col) = 0 ;
387  for (int col=nr; col<nrtot; col++)
388  oper.set(nr+nr_eq2, col) = 1 ;
389  for (int col=0; col<nr; col++)
390  oper.set(nr+nr_eq2+1, col) = 0 ;
391  pari = -1 ;
392  for (int col=0; col<nr; col++) {
393  oper.set(nr+nr_eq2+1, nr+col) = pari*col*col ;
394  pari = pari ;
395  }
396  if (l_q > 1) {
397  for (int col=0; col<nr; col++)
398  oper.set(nr+nr_eq2+2, col) = 0 ;
399  oper.set(nr+nr_eq2+2, nr) = 1 ;
400  for (int col=nr+1; col<nrtot; col++)
401  oper.set(nr+nr_eq2+2, col) = 0 ;
402  }
403  oper.set_lu() ;
404 
405  // Filling the r.h.s
406  //------------------
407  for (int i=0; i<nr_eq1; i++)
408  sec_membre.set(i) = (*S_eta.get_spectral_va().c_cf)(nzm1, k, j, i) ;
409  for (int i=nr_eq1; i<nr; i++)
410  sec_membre.set(i) = 0 ;
411  for (int i=0; i<nr_eq2; i++)
412  sec_membre.set(i+nr) =(*S_r.get_spectral_va().c_cf)(nzm1, k, j, i);
413  for (int i=nr_eq2; i<nr; i++)
414  sec_membre.set(nr+i) = 0 ;
415  Tbl big_res = oper.inverse(sec_membre) ;
416 
417  // Putting coefficients of h and v to individual arrays
418  //-----------------------------------------------------
419  for (int i=0; i<nr; i++) {
420  sol_part_eta.set(nzm1, k, j, i) = big_res(i) ;
421  sol_part_vr.set(nzm1, k, j, i) = big_res(i+nr) ;
422  }
423 
424  // Homogeneous solutions (only 1/r^(l+2) and 1/r^l in the CED)
425  //------------------------------------------------------------
426  if (l_q == 1) {
427  Tbl sol_hom1 = solh(nr, 0, 0., base_r) ;
428  Tbl sol_hom2 = solh(nr, 2, 0., base_r) ;
429  double fac_eta = lam + 2. ;
430  double fac_vr = 2*lam + 2. ;
431  for (int i=0 ; i<nr ; i++) {
432  sol_hom_deux_eta.set(nzm1, k, j, i) = sol_hom2(i) ;
433  sol_hom_quatre_eta.set(nzm1, k, j, i) = fac_eta*sol_hom1(i) ;
434  sol_hom_deux_vr.set(nzm1, k, j, i) = -2*sol_hom2(i) ;
435  sol_hom_quatre_vr.set(nzm1, k, j, i) = fac_vr*sol_hom1(i) ;
436  }
437  }
438  else {
439  sec_membre.annule_hard() ;
440  sec_membre.set(nr-1) = 1 ;
441  big_res = oper.inverse(sec_membre) ;
442 
443  for (int i=0; i<nr; i++) {
444  sol_hom_deux_eta.set(nzm1, k, j, i) = big_res(i) ;
445  sol_hom_deux_vr.set(nzm1, k, j, i) = big_res(nr+i) ;
446  }
447  sec_membre.set(nr-1) = 0 ;
448  sec_membre.set(2*nr-1) = 1 ;
449  big_res = oper.inverse(sec_membre) ;
450  for (int i=0; i<nr; i++) {
451  sol_hom_quatre_eta.set(nzm1, k, j, i) = big_res(i) ;
452  sol_hom_quatre_vr.set(nzm1, k, j, i) = big_res(nr+i) ;
453  }
454  }
455  }
456  }
457  }
458 
459  // Now let's match everything ...
460  //-------------------------------
461 
462  // Resulting V^r & eta
463  vr.set_etat_qcq() ;
464  vr.set_spectral_base(base) ;
466  Mtbl_cf& cf_vr = *vr.set_spectral_va().c_cf ;
467  cf_vr.annule_hard() ;
468  het.set_etat_qcq() ;
469  het.set_spectral_base(base) ;
471  Mtbl_cf& cf_eta = *het.set_spectral_va().c_cf ;
472  cf_eta.annule_hard() ;
473  int taille = 4*nzm1 -2 ; // Two less than R3 case
474  Tbl sec_membre(taille) ;
475  Matrice systeme(taille, taille) ;
476  systeme.set_etat_qcq() ;
477  int ligne ; int colonne ;
478 
479  // Loop on l and m
480  //----------------
481  for (int k=0 ; k<np+1 ; k++)
482  for (int j=0 ; j<nt ; j++) {
483  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
484  if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q != 0)) {
485 
486  ligne = 0 ;
487  colonne = 0 ;
488  systeme.annule_hard() ;
489  sec_membre.annule_hard() ;
490 
491 
492  //shell 1 (boundary condition)
493  int zone1 = 1;
494  nr = mg.get_nr(zone1) ;
495  alpha = mpaff->get_alpha()[zone1] ;
496 
497  // Here we prepare for a Robyn-type boundary condition on eta.
498  // Parameters are dir_eta and neum_eta
499  systeme.set(ligne, colonne)
500  = -dir_eta*sol_hom_un_eta.val_in_bound_jk(zone1, j, k) ;
501  systeme.set(ligne, colonne+1)
502  = -dir_eta * sol_hom_deux_eta.val_in_bound_jk(zone1, j, k) ;
503  systeme.set(ligne, colonne+2)
504  = -dir_eta*sol_hom_trois_eta.val_in_bound_jk(zone1, j, k) ;
505  systeme.set(ligne, colonne+3)
506  = -dir_eta*sol_hom_quatre_eta.val_in_bound_jk(zone1, j, k) ;
507 
508 
509 
510 
511  sec_membre.set(ligne) += dir_eta* sol_part_eta.val_in_bound_jk(zone1, j, k) ;
512 
513 
514 
515  int pari = -1 ;
516  for (int i=0; i<nr; i++) {
517  systeme.set(ligne, colonne)
518  -= neum_eta*pari*i*i*sol_hom_un_eta(zone1, k, j, i)/alpha ;
519  systeme.set(ligne, colonne+1)
520  -= neum_eta*pari*i*i*sol_hom_deux_eta(zone1, k, j, i)/alpha ;
521  systeme.set(ligne, colonne+2)
522  -= neum_eta*pari*i*i*sol_hom_trois_eta(zone1, k, j, i)/alpha ;
523  systeme.set(ligne, colonne+3)
524  -= neum_eta*pari*i*i*sol_hom_quatre_eta(zone1, k, j, i)/alpha ;
525  sec_membre.set(ligne)
526  += neum_eta*pari*i*i* sol_part_eta(zone1, k, j, i)/alpha ;
527  pari = -pari ;
528  }
529 
530  sec_membre.set(ligne) -= (*bound_eta).val_in_bound_jk(1,j,k);
531  ligne++ ;
532 
533 
534 
535  // ... and their couterparts for V^r
536  systeme.set(ligne, colonne)
537  = - dir_vr*sol_hom_un_vr.val_in_bound_jk(zone1, j, k) ;
538  systeme.set(ligne, colonne+1)
539  = - dir_vr*sol_hom_deux_vr.val_in_bound_jk(zone1, j, k) ;
540  systeme.set(ligne, colonne+2)
541  = -dir_vr*sol_hom_trois_vr.val_in_bound_jk(zone1, j, k) ;
542  systeme.set(ligne, colonne+3)
543  = -dir_vr*sol_hom_quatre_vr.val_in_bound_jk(zone1, j, k) ;
544  sec_membre.set(ligne) += dir_vr*sol_part_vr.val_in_bound_jk(zone1, j, k) ;
545 
546 
547 
548 
549  pari = -1 ;
550  for (int i=0; i<nr; i++) {
551  systeme.set(ligne, colonne)
552  -= neum_vr*pari*i*i*sol_hom_un_vr(zone1, k, j, i)/alpha ;
553  systeme.set(ligne, colonne+1)
554  -= neum_vr*pari*i*i*sol_hom_deux_vr(zone1, k, j, i)/alpha ;
555  systeme.set(ligne, colonne+2)
556  -= neum_vr*pari*i*i*sol_hom_trois_vr(zone1, k, j, i)/alpha ;
557  systeme.set(ligne, colonne+3)
558  -= neum_vr*pari*i*i*sol_hom_quatre_vr(zone1, k, j, i)/alpha ;
559  sec_membre.set(ligne)
560  += neum_vr*pari*i*i* sol_part_vr(zone1, k, j, i)/alpha ;
561  pari = -pari ;
562  }
563 
564  sec_membre.set(ligne) -= (*bound_vr).val_in_bound_jk(1,j,k);
565  ligne++ ;
566 
567 
568 
569 
570  // Values in 1: beginning with solution in eta...
571 
572  systeme.set(ligne, colonne)
573  += sol_hom_un_eta.val_out_bound_jk(zone1, j, k) ;
574  systeme.set(ligne, colonne+1)
575  += sol_hom_deux_eta.val_out_bound_jk(zone1, j, k) ;
576  systeme.set(ligne, colonne+2)
577  += sol_hom_trois_eta.val_out_bound_jk(zone1, j, k) ;
578  systeme.set(ligne, colonne+3)
579  += sol_hom_quatre_eta.val_out_bound_jk(zone1, j, k) ;
580  sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone1, j, k) ;
581  ligne++ ;
582  // ... and their counterparts for V^r
583  systeme.set(ligne, colonne)
584  += sol_hom_un_vr.val_out_bound_jk(zone1, j, k) ;
585  systeme.set(ligne, colonne+1)
586  += sol_hom_deux_vr.val_out_bound_jk(zone1, j, k) ;
587  systeme.set(ligne, colonne+2)
588  += sol_hom_trois_vr.val_out_bound_jk(zone1, j, k) ;
589  systeme.set(ligne, colonne+3)
590  += sol_hom_quatre_vr.val_out_bound_jk(zone1, j, k) ;
591  sec_membre.set(ligne) -= sol_part_vr.val_out_bound_jk(zone1, j, k) ;
592  ligne++ ;
593 
594  //derivative at 1
595  pari = 1 ;
596  for (int i=0; i<nr; i++) {
597  systeme.set(ligne, colonne)
598  += pari*i*i*sol_hom_un_eta(zone1, k, j, i)/alpha ;
599  systeme.set(ligne, colonne+1)
600  += pari*i*i*sol_hom_deux_eta(zone1, k, j, i)/alpha ;
601  systeme.set(ligne, colonne+2)
602  += pari*i*i*sol_hom_trois_eta(zone1, k, j, i)/alpha ;
603  systeme.set(ligne, colonne+3)
604  += pari*i*i*sol_hom_quatre_eta(zone1, k, j, i)/alpha ;
605  sec_membre.set(ligne)
606  -= pari*i*i* sol_part_eta(zone1, k, j, i)/alpha ;
607  pari = pari ;
608  }
609  ligne++ ;
610  // ... and their counterparts for V^r
611  pari = 1 ;
612  for (int i=0; i<nr; i++) {
613  systeme.set(ligne, colonne)
614  += pari*i*i*sol_hom_un_vr(zone1, k, j, i)/alpha ;
615  systeme.set(ligne, colonne+1)
616  += pari*i*i*sol_hom_deux_vr(zone1, k, j, i)/alpha ;
617  systeme.set(ligne, colonne+2)
618  += pari*i*i*sol_hom_trois_vr(zone1, k, j, i)/alpha ;
619  systeme.set(ligne, colonne+3)
620  += pari*i*i*sol_hom_quatre_vr(zone1, k, j, i)/alpha ;
621  sec_membre.set(ligne)
622  -= pari*i*i* sol_part_vr(zone1, k, j, i)/alpha ;
623  pari = pari ;
624  }
625  colonne += 4 ;
626 
627  // Other shells
628  if ( nz <= 2){
629  cout <<"WARNING!! Mapping must contain at least 2 shells!!" << endl;}
630  for (int zone=2 ; zone<nzm1 ; zone++) {
631 
632  nr = mg.get_nr(zone) ;
633  alpha = mpaff->get_alpha()[zone] ;
634  ligne -= 3 ;
635  systeme.set(ligne, colonne)
636  = -sol_hom_un_eta.val_in_bound_jk(zone, j, k) ;
637  systeme.set(ligne, colonne+1)
638  = -sol_hom_deux_eta.val_in_bound_jk(zone, j, k) ;
639  systeme.set(ligne, colonne+2)
640  = -sol_hom_trois_eta.val_in_bound_jk(zone, j, k) ;
641  systeme.set(ligne, colonne+3)
642  = -sol_hom_quatre_eta.val_in_bound_jk(zone, j, k) ;
643  sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(zone, j, k) ;
644  ligne++ ;
645  // ... and their counterparts for V^r
646  systeme.set(ligne, colonne)
647  = -sol_hom_un_vr.val_in_bound_jk(zone, j, k) ;
648  systeme.set(ligne, colonne+1)
649  = -sol_hom_deux_vr.val_in_bound_jk(zone, j, k) ;
650  systeme.set(ligne, colonne+2)
651  = -sol_hom_trois_vr.val_in_bound_jk(zone, j, k) ;
652  systeme.set(ligne, colonne+3)
653  = -sol_hom_quatre_vr.val_in_bound_jk(zone, j, k) ;
654  sec_membre.set(ligne) += sol_part_vr.val_in_bound_jk(zone, j, k) ;
655  ligne++ ;
656 
657  //derivative of (x+echelle)^(l-1) at -1
658  pari = -1 ;
659  for (int i=0; i<nr; i++) {
660  systeme.set(ligne, colonne)
661  -= pari*i*i*sol_hom_un_eta(zone, k, j, i)/alpha ;
662  systeme.set(ligne, colonne+1)
663  -= pari*i*i*sol_hom_deux_eta(zone, k, j, i)/alpha ;
664  systeme.set(ligne, colonne+2)
665  -= pari*i*i*sol_hom_trois_eta(zone, k, j, i)/alpha ;
666  systeme.set(ligne, colonne+3)
667  -= pari*i*i*sol_hom_quatre_eta(zone, k, j, i)/alpha ;
668  sec_membre.set(ligne)
669  += pari*i*i* sol_part_eta(zone, k, j, i)/alpha ;
670  pari = -pari ;
671  }
672  ligne++ ;
673  // ... and their counterparts for V^r
674  pari = -1 ;
675  for (int i=0; i<nr; i++) {
676  systeme.set(ligne, colonne)
677  -= pari*i*i*sol_hom_un_vr(zone, k, j, i)/alpha ;
678  systeme.set(ligne, colonne+1)
679  -= pari*i*i*sol_hom_deux_vr(zone, k, j, i)/alpha ;
680  systeme.set(ligne, colonne+2)
681  -= pari*i*i*sol_hom_trois_vr(zone, k, j, i)/alpha ;
682  systeme.set(ligne, colonne+3)
683  -= pari*i*i*sol_hom_quatre_vr(zone, k, j, i)/alpha ;
684  sec_membre.set(ligne)
685  += pari*i*i* sol_part_vr(zone, k, j, i)/alpha ;
686  pari = -pari ;
687  }
688  ligne++ ;
689 
690  systeme.set(ligne, colonne)
691  += sol_hom_un_eta.val_out_bound_jk(zone, j, k) ;
692  systeme.set(ligne, colonne+1)
693  += sol_hom_deux_eta.val_out_bound_jk(zone, j, k) ;
694  systeme.set(ligne, colonne+2)
695  += sol_hom_trois_eta.val_out_bound_jk(zone, j, k) ;
696  systeme.set(ligne, colonne+3)
697  += sol_hom_quatre_eta.val_out_bound_jk(zone, j, k) ;
698  sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone, j, k) ;
699  ligne++ ;
700  // ... and their counterparts for V^r
701  systeme.set(ligne, colonne)
702  += sol_hom_un_vr.val_out_bound_jk(zone, j, k) ;
703  systeme.set(ligne, colonne+1)
704  += sol_hom_deux_vr.val_out_bound_jk(zone, j, k) ;
705  systeme.set(ligne, colonne+2)
706  += sol_hom_trois_vr.val_out_bound_jk(zone, j, k) ;
707  systeme.set(ligne, colonne+3)
708  += sol_hom_quatre_vr.val_out_bound_jk(zone, j, k) ;
709  sec_membre.set(ligne) -= sol_part_vr.val_out_bound_jk(zone, j, k) ;
710  ligne++ ;
711 
712  //derivative at 1
713  pari = 1 ;
714  for (int i=0; i<nr; i++) {
715  systeme.set(ligne, colonne)
716  += pari*i*i*sol_hom_un_eta(zone, k, j, i)/alpha ;
717  systeme.set(ligne, colonne+1)
718  += pari*i*i*sol_hom_deux_eta(zone, k, j, i)/alpha ;
719  systeme.set(ligne, colonne+2)
720  += pari*i*i*sol_hom_trois_eta(zone, k, j, i)/alpha ;
721  systeme.set(ligne, colonne+3)
722  += pari*i*i*sol_hom_quatre_eta(zone, k, j, i)/alpha ;
723  sec_membre.set(ligne)
724  -= pari*i*i* sol_part_eta(zone, k, j, i)/alpha ;
725  pari = pari ;
726  }
727  ligne++ ;
728  // ... and their couterparts for V^r
729  pari = 1 ;
730  for (int i=0; i<nr; i++) {
731  systeme.set(ligne, colonne)
732  += pari*i*i*sol_hom_un_vr(zone, k, j, i)/alpha ;
733  systeme.set(ligne, colonne+1)
734  += pari*i*i*sol_hom_deux_vr(zone, k, j, i)/alpha ;
735  systeme.set(ligne, colonne+2)
736  += pari*i*i*sol_hom_trois_vr(zone, k, j, i)/alpha ;
737  systeme.set(ligne, colonne+3)
738  += pari*i*i*sol_hom_quatre_vr(zone, k, j, i)/alpha ;
739  sec_membre.set(ligne)
740  -= pari*i*i* sol_part_vr(zone, k, j, i)/alpha ;
741  pari = pari ;
742  }
743  colonne += 4 ;
744  }
745 
746  //Compactified external domain
747  nr = mg.get_nr(nzm1) ;
748 
749  alpha = mpaff->get_alpha()[nzm1] ;
750  ligne -= 3 ;
751  //value of (x-1)^(l+2) at -1 :
752  systeme.set(ligne, colonne)
753  -= sol_hom_deux_eta.val_in_bound_jk(nzm1, j, k) ;
754  systeme.set(ligne, colonne+1)
755  -= sol_hom_quatre_eta.val_in_bound_jk(nzm1, j, k) ;
756  sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(nzm1, j, k) ;
757  //... and of its counterpart for V^r
758  systeme.set(ligne+1, colonne)
759  -= sol_hom_deux_vr.val_in_bound_jk(nzm1, j, k) ;
760  systeme.set(ligne+1, colonne+1)
761  -= sol_hom_quatre_vr.val_in_bound_jk(nzm1, j, k) ;
762  sec_membre.set(ligne+1) += sol_part_vr.val_in_bound_jk(nzm1, j, k) ;
763 
764  ligne += 2 ;
765  //derivative of (x-1)^(l+2) at -1 :
766  pari = 1 ;
767  for (int i=0; i<nr; i++) {
768  systeme.set(ligne, colonne)
769  -= 4*alpha*pari*i*i*sol_hom_deux_eta(nzm1, k, j, i) ;
770  systeme.set(ligne, colonne+1)
771  -= 4*alpha*pari*i*i*sol_hom_quatre_eta(nzm1, k, j, i) ;
772  sec_membre.set(ligne)
773  += 4*alpha*pari*i*i* sol_part_eta(nzm1, k, j, i) ;
774  pari = -pari ;
775  }
776  ligne++ ;
777  // ... and their counterparts for V^r
778  pari = 1 ;
779  for (int i=0; i<nr; i++) {
780  systeme.set(ligne, colonne)
781  -= 4*alpha*pari*i*i*sol_hom_deux_vr(nzm1, k, j, i) ;
782  systeme.set(ligne, colonne+1)
783  -= 4*alpha*pari*i*i*sol_hom_quatre_vr(nzm1, k, j, i) ;
784  sec_membre.set(ligne)
785  += 4*alpha*pari*i*i* sol_part_vr(nzm1, k, j, i) ;
786  pari = -pari ;
787  }
788 
789  // Solution of the system giving the coefficients for the homogeneous
790  // solutions
791  //-------------------------------------------------------------------
792  systeme.set_lu() ;
793  Tbl facteurs(systeme.inverse(sec_membre)) ;
794  int conte = 0 ;
795 
796  // everything is put to the right place, the same combination of hom.
797  // solutions (with some l or -(l+1) factors) must be used for V^r
798  //-------------------------------------------------------------------
799  nr = mg.get_nr(0) ; //nucleus
800  for (int i=0 ; i<nr ; i++) {
801  cf_eta.set(0, k, j, i) = 0.;
802  cf_vr.set(0, k, j, i) = 0.;
803  }
804 
805  for (int zone=1 ; zone<nzm1 ; zone++) { //shells
806  nr = mg.get_nr(zone) ;
807  for (int i=0 ; i<nr ; i++) {
808  cf_eta.set(zone, k, j, i) =
809  sol_part_eta(zone, k, j, i)
810  +facteurs(conte)*sol_hom_un_eta(zone, k, j, i)
811  +facteurs(conte+1)*sol_hom_deux_eta(zone, k, j, i)
812  +facteurs(conte+2)*sol_hom_trois_eta(zone, k, j, i)
813  +facteurs(conte+3)*sol_hom_quatre_eta(zone, k, j, i) ;
814  cf_vr.set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
815  +facteurs(conte)*sol_hom_un_vr(zone, k, j, i)
816  +facteurs(conte+1)*sol_hom_deux_vr(zone, k, j, i)
817  +facteurs(conte+2)*sol_hom_trois_vr(zone, k, j, i)
818  +facteurs(conte+3)*sol_hom_quatre_vr(zone, k, j, i) ;
819  }
820  conte+=4 ;
821  }
822  nr = mg.get_nr(nz-1) ; //compactified external domain
823  for (int i=0 ; i<nr ; i++) {
824  cf_eta.set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
825  +facteurs(conte)*sol_hom_deux_eta(nzm1, k, j, i)
826  +facteurs(conte+1)*sol_hom_quatre_eta(nzm1, k, j, i) ;
827  cf_vr.set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
828  +facteurs(conte)*sol_hom_deux_vr(nzm1, k, j, i)
829  +facteurs(conte+1)*sol_hom_quatre_vr(nzm1, k, j, i) ;
830 
831  }
832  } // End of nullite_plm
833  } //End of loop on theta
834  vr.set_spectral_va().ylm_i() ;
835  vr += vrl0 ;
836  het.set_spectral_va().ylm_i() ;
837 
838 
839 
840  Scalar pipo(*mpaff);
841  pipo.annule_hard();
842  pipo.std_spectral_base();
843  pipo.set_dzpuis(3);
844  Param_elliptic param_mu(mu());
845 
846 
847 
848 
849  resu.set_vr_eta_mu(vr, het, mu().sol_elliptic_boundary(param_mu, (*bound_mu), dir_mu , neum_mu)) ;
850 
851  return ;
852 
853 }
854 }
Lorene::Diff_sxdsdx::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_sxdsdx.C:97
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Lorene::Base_val
Bases of the spectral expansions.
Definition: base_val.h:322
Lorene::Diff_sx2
Class for the elementary differential operator division by (see the base class Diff ).
Definition: diff.h:369
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Lorene::Valeur::set_etat_cf_qcq
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
Definition: valeur.C:712
Lorene::Matrice::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:175
Lorene::Diff_sxdsdx
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:450
Lorene::Scalar::set_spectral_base
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
Definition: scalar.C:797
Lorene::Diff_xdsdx2::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_xdsdx2.C:97
Lorene::Mg3d
Multi-domain grid.
Definition: grilles.h:273
Lorene::Vector::poisson_boundary2
void poisson_boundary2(double lam, Vector &resu, Scalar boundvr, Scalar boundeta, Scalar boundmu, double dir_vr, double neum_vr, double dir_eta, double neum_eta, double dir_mu, double neum_mu) const
Alternative to previous poisson_boundary method for vectors ; this uses method 6 for vectorial solvin...
Definition: vector_poisson_boundary2.C:80
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:372
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Tensor::cmp
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition: tensor.h:315
Lorene::Mg3d::get_type_r
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:474
Lorene::Diff_dsdx2::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_dsdx2.C:91
Lorene::Matrice::annule_hard
void annule_hard()
Sets the logical state to ETATQCQ (undefined state).
Definition: matrice.C:193
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Mtbl_cf::set
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:294
Lorene::Param_elliptic
This class contains the parameters needed to call the general elliptic solver.
Definition: param_elliptic.h:135
Lorene::Scalar::get_spectral_base
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va
Definition: scalar.h:1294
Lorene::Tbl::t
double * t
The array of double.
Definition: tbl.h:173
Lorene::Scalar::get_spectral_va
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:601
Lorene::Matrice::inverse
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition: matrice.C:424
Lorene::Mtbl_cf::val_in_bound_jk
double val_in_bound_jk(int l, int j, int k) const
Computes the angular coefficient of index j,k of the field represented by *this at by means of the s...
Definition: mtbl_cf_val_point.C:452
Lorene::Diff_xdsdx::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_xdsdx.C:98
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
Lorene::Base_val::give_quant_numbers
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Definition: base_val_quantum.C:81
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Lorene::Diff_xdsdx2
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:531
Lorene::Vector::set_vr_eta_mu
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
Lorene::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:295
Lorene::Base_vect_spher
Spherical orthonormal vectorial bases (triads).
Definition: base_vect.h:308
Lorene::Mtbl_cf::annule_hard
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition: mtbl_cf.C:312
Lorene::Vector::mu
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written:
Definition: vector_etamu.C:98
Lorene::Diff_sx2::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_sx2.C:98
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Diff_xdsdx
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:409
Lorene::Tensor::triad
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tensor.h:303
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Lorene::Map_af
Affine radial mapping.
Definition: map.h:2027
Lorene::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:353
Lorene::Diff_id
Class for the elementary differential operator Identity (see the base class Diff ).
Definition: diff.h:210
Lorene::Mtbl_cf::val_out_bound_jk
double val_out_bound_jk(int l, int j, int k) const
Computes the angular coefficient of index j,k of the field represented by *this at by means of the s...
Definition: mtbl_cf_val_point.C:428
Lorene::Diff_x2dsdx2
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:490
Lorene::Scalar::set_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:808
Lorene::Vector
Tensor field of valence 1.
Definition: vector.h:188
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
Lorene::Vector::eta
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written:
Definition: vector_etamu.C:66
Lorene::Param_elliptic::set_poisson_vect_r
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.
Definition: param_elliptic.C:487
Lorene::Matrice
Matrix handling.
Definition: matrice.h:152
Lorene::Map_af::get_beta
const double * get_beta() const
Returns the pointer on the array beta.
Definition: map_af.C:481
Lorene::Scalar::std_spectral_base
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition: scalar.C:784
Lorene::Scalar::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:380
Lorene::Diff_dsdx2
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:172
Lorene::Valeur::base
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
Lorene::Diff_dsdx::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_dsdx.C:94
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
Lorene::Matrice::set
double & set(int j, int i)
Read/write of a particuliar element.
Definition: matrice.h:277
Lorene::Diff_x2dsdx2::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_x2dsdx2.C:95
Lorene::Diff_dsdx
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:129
Lorene::Map_af::get_alpha
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:477
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:186
Lorene::Scalar::sol_elliptic_boundary
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
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131
Lorene::Matrice::set_lu
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition: matrice.C:392