LORENE
tensorellipticBt.C
1 #include<math.h>
2 // Lorene headers
3 #include "metric.h"
4 #include "nbr_spx.h"
5 #include "utilitaires.h"
6 #include "graphique.h"
7 #include "proto.h"
8 #include "diff.h"
9 
10 
11 namespace Lorene {
12 // Inversion of the weakly degenerate eliptic operator associatd with spectral quantity tilde(B), with main characteristics
13 //fit and fit2 (Suitable for tensorial resolution of BH spacetime)
14 //See also function tilde(laplacian)
15 
16  void tensorellipticBt( Scalar source, Scalar& resu, double fit, double fit2, double fit0d2, double fit1d2, double fit0d3, double fit1d3) {
17 
18 
19  const int nz = (*source.get_mp().get_mg()).get_nzone(); // Number of domains
20  int nr = (*source.get_mp().get_mg()).get_nr(1); // Number of collocation points in r in each domain
21  int nt = (*source.get_mp().get_mg()).get_nt(1); // Number of collocation points in theta in each domain
22  int np = (*source.get_mp().get_mg()).get_np(1); // Number of collocation points in phi in each domain
23  const Map_af* map = dynamic_cast<const Map_af*>(&source.get_mp()) ;
24  const Mg3d* mgrid = (*map).get_mg();
25 
26 
27  // Some helpful stuff...
28 
29  const Coord& rr = (*map).r ;
30  Scalar rrr (*map) ;
31  rrr = rr ;
32  rrr.set_spectral_va().set_base(source.get_spectral_va().base);
33 
34  Scalar source_coq = source ;
35  source_coq.mult_r() ;
36  source_coq.mult_r() ;
37  source.set_spectral_va().ylm() ;
38  source_coq.set_spectral_va().ylm() ;
39  Scalar phi(source.get_mp()) ;
40  phi.annule_hard() ;
41  // phi.std_spectral_base();
42  phi.set_spectral_va().set_base(source.get_spectral_va().base) ;
43  phi.set_spectral_va().ylm() ;
44  Mtbl_cf& sol_coef = (*phi.set_spectral_va().c_cf) ;
45 
46  const Base_val& base = source.get_spectral_base() ;
47  Mtbl_cf sol_part(mgrid, base) ; sol_part.annule_hard() ;
48  Mtbl_cf sol_hom1(mgrid, base) ; sol_hom1.annule_hard() ;
49  Mtbl_cf sol_hom2(mgrid, base) ; sol_hom2.annule_hard() ;
50 
51  int l_q, m_q, base_r ;
52 
53  { int lz = 0 ;
54  for (int k=0; k < np; k++)
55  for (int j=0; j<nt; j++) {
56  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
57  if (nullite_plm(j, nt, k, np, base) == 1) {
58  for (int ii=0 ; ii<nr ; ii++){
59  sol_hom1.set(lz, k, j, ii) = 0 ;
60  sol_part.set(lz, k, j, ii) = 0 ;
61  }
62 
63  }
64  }
65  }
66 
67 
68  { int lz = 1 ; // The first shell is a really particular case, where the operator is different, and homogeneous solutions have to be handled really carefully.
69  double alpha = (*map).get_alpha()[lz] ;
70  double beta = (*map).get_beta()[lz] ;
71  double ech = beta / alpha ;
72  for (int k=0; k < np; k++)
73  for (int j=0; j<nt; j++) {
74  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
75  if (nullite_plm(j, nt, k, np, base) == 1) {
76 
77 
78 
79  Matrice ope(nr,nr) ;
80  ope.annule_hard() ;
81 
82  Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
83  Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
84  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
85  Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
86  Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
87  Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
88  Diff_x3dsdx2 x3dx2 (base_r, nr); const Matrice& mx3dx2 = x3dx2.get_matrice();
89  Diff_x4dsdx2 x4dx2 (base_r, nr); const Matrice& mx4dx2 = x4dx2.get_matrice();
90 
91 
92 
93 // ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
94 // - l_q*(l_q+1)*mid - (mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) - (fit)*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) + (fit)*beta*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) + (fit)*alpha*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2) ;
95 
96 
97 // ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
98 // - l_q*(l_q+1)*mid - ((mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) -(fit)*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) + (fit)*beta*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2)
99 // + (fit)*alpha*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2) + fit2*alpha*alpha*(mx4dx2 + 2*ech*mx3dx2 + ech*ech*mx2dx2)
100 // + 2*fit2*alpha*(beta-1.)*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2) + fit2*(beta-1.)*(beta-1.)*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2));
101 
102 
103 
104  ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
105  - l_q*(l_q+1)*mid + 2*l_q*mid - ((mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) -(fit)*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) + (fit)*beta*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2)
106  + (fit)*alpha*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2) + fit2*alpha*alpha*(mx4dx2 + 2*ech*mx3dx2 + ech*ech*mx2dx2)
107  + fit2*alpha*(beta-1.)*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2) + fit2*alpha*(beta- rrr.val_grid_point(1,0,0, nr-1))*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2)+ fit2*(beta-1.)*(beta- rrr.val_grid_point(1,0,0,nr -1))*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2));
108 
109 
110 
111 
112  for (int col=0; col<nr; col++)
113  ope.set(nr-1, col) = 0 ;
114  ope.set(nr-1, 0) = 1 ;
115 
116 
117  Tbl rhs(nr);
118  rhs.annule_hard() ;
119  for (int i=0; i<nr; i++)
120  rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(1, k, j, i) ;
121  rhs.set(nr-1) = 0 ;
122  ope.set_lu() ;
123  Tbl sol = ope.inverse(rhs) ;
124 
125 
126  for (int i=0; i<nr; i++)
127  sol_part.set(1, k, j, i) = sol(i) ;
128 
129  rhs.annule_hard();
130  rhs.set(nr-1) = 1 ;
131  sol = ope.inverse(rhs) ;
132 
133 
134  for (int i=0; i<nr; i++)
135  sol_hom1.set(1, k, j, i) = sol(i) ;
136 
137 
138  }
139  }
140  }
141 
142  // Attention! zones 2 et 3traitee separement egalement!!!
143 
144 
145  // Current implementations only allow grids with more than 3 shells.
146 
147  { int lz = 2 ;
148 
149  double alpha = (*map).get_alpha()[lz] ;
150  double beta = (*map).get_beta()[lz] ;
151  double ech = beta / alpha ;
152 
153  for (int k=0; k < np; k++)
154  for (int j=0; j<nt; j++) {
155  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
156  if (nullite_plm(j, nt, k, np, base) == 1) {
157 
158  Matrice ope(nr,nr) ;
159  ope.annule_hard() ;
160 
161  Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
162  Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
163  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
164  Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
165  Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
166  Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
167  Diff_x3dsdx2 x3dx2 (base_r, nr); const Matrice& mx3dx2 = x3dx2.get_matrice();
168 
169  ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
170  - l_q*(l_q+1)*mid + 2*l_q*mid - fit0d2*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) + (fit1d2)*(rrr.val_grid_point(lz, 0, 0, 0))*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) - (fit1d2)*beta*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) - (fit1d2)*alpha*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2) ;
171 
172 
173 
174 // ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
175 // - l_q*(l_q+1)*mid ;
176 
177  for (int col=0; col<nr; col++)
178  ope.set(nr-1, col) = 0 ;
179  ope.set(nr-1, 0) = 1 ;
180  for (int col=0; col<nr; col++) {
181  ope.set(nr-2, col) = 0 ;
182  }
183  ope.set(nr-2, 1) = 1 ;
184 
185  Tbl rhs(nr) ;
186  rhs.annule_hard() ;
187  for (int i=0; i<nr; i++)
188  rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(lz, k, j, i) ;
189  rhs.set(nr-2) = 0 ;
190  rhs.set(nr-1) = 0 ;
191  ope.set_lu() ;
192  Tbl sol = ope.inverse(rhs) ;
193 
194  for (int i=0; i<nr; i++)
195  sol_part.set(lz, k, j, i) = sol(i) ;
196 
197  rhs.annule_hard() ;
198  rhs.set(nr-2) = 1 ;
199  sol = ope.inverse(rhs) ;
200  for (int i=0; i<nr; i++)
201  sol_hom1.set(lz, k, j, i) = sol(i) ;
202 
203  rhs.set(nr-2) = 0 ;
204  rhs.set(nr-1) = 1 ;
205  sol = ope.inverse(rhs) ;
206  for (int i=0; i<nr; i++)
207  sol_hom2.set(lz, k, j, i) = sol(i) ;
208 
209  }
210  }
211 
212  }
213 
214 
215 
216  { int lz = 3 ;
217 
218  double alpha = (*map).get_alpha()[lz] ;
219  double beta = (*map).get_beta()[lz] ;
220  double ech = beta / alpha ;
221 
222  for (int k=0; k < np; k++)
223  for (int j=0; j<nt; j++) {
224  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
225  if (nullite_plm(j, nt, k, np, base) == 1) {
226 
227  Matrice ope(nr,nr) ;
228  ope.annule_hard() ;
229 
230  Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
231  Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
232  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
233  Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
234  Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
235  Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
236  Diff_x3dsdx2 x3dx2 (base_r, nr); const Matrice& mx3dx2 = x3dx2.get_matrice();
237 
238  ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
239  - l_q*(l_q+1)*mid +2*l_q*mid - fit0d3*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) + (fit1d3)*(rrr.val_grid_point(lz, 0, 0, 0))*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) - (fit1d3)*beta*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2) - (fit1d3)*alpha*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2) ;
240 
241 
242 
243 // ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
244 // - l_q*(l_q+1)*mid ;
245 
246  for (int col=0; col<nr; col++)
247  ope.set(nr-1, col) = 0 ;
248  ope.set(nr-1, 0) = 1 ;
249  for (int col=0; col<nr; col++) {
250  ope.set(nr-2, col) = 0 ;
251  }
252  ope.set(nr-2, 1) = 1 ;
253 
254  Tbl rhs(nr) ;
255  rhs.annule_hard() ;
256  for (int i=0; i<nr; i++)
257  rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(lz, k, j, i) ;
258  rhs.set(nr-2) = 0 ;
259  rhs.set(nr-1) = 0 ;
260  ope.set_lu() ;
261  Tbl sol = ope.inverse(rhs) ;
262 
263  for (int i=0; i<nr; i++)
264  sol_part.set(lz, k, j, i) = sol(i) ;
265 
266  rhs.annule_hard() ;
267  rhs.set(nr-2) = 1 ;
268  sol = ope.inverse(rhs) ;
269  for (int i=0; i<nr; i++)
270  sol_hom1.set(lz, k, j, i) = sol(i) ;
271 
272  rhs.set(nr-2) = 0 ;
273  rhs.set(nr-1) = 1 ;
274  sol = ope.inverse(rhs) ;
275  for (int i=0; i<nr; i++)
276  sol_hom2.set(lz, k, j, i) = sol(i) ;
277 
278  }
279  }
280 
281  }
282 
283 
284 
285 
286 
287 
288 
289  // Current implementations only allow grids with more than 2 shells.
290 
291  for (int lz=4; lz<nz-1; lz++) {
292  double ech = (*map).get_beta()[lz] / (*map).get_alpha()[lz] ;
293  for (int k=0; k < np; k++)
294  for (int j=0; j<nt; j++) {
295  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
296  if (nullite_plm(j, nt, k, np, base) == 1) {
297 
298  Matrice ope(nr,nr) ;
299  ope.annule_hard() ;
300 
301  Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
302  Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
303  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
304  Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
305  Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
306  Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
307  ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
308  - l_q*(l_q+1)*mid +2*l_q*mid ;
309 
310  for (int col=0; col<nr; col++)
311  ope.set(nr-1, col) = 0 ;
312  ope.set(nr-1, 0) = 1 ;
313  for (int col=0; col<nr; col++) {
314  ope.set(nr-2, col) = 0 ;
315  }
316  ope.set(nr-2, 1) = 1 ;
317 
318  Tbl rhs(nr) ;
319  rhs.annule_hard() ;
320  for (int i=0; i<nr; i++)
321  rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(lz, k, j, i) ;
322  rhs.set(nr-2) = 0 ;
323  rhs.set(nr-1) = 0 ;
324  ope.set_lu() ;
325  Tbl sol = ope.inverse(rhs) ;
326 
327  for (int i=0; i<nr; i++)
328  sol_part.set(lz, k, j, i) = sol(i) ;
329 
330  rhs.annule_hard() ;
331  rhs.set(nr-2) = 1 ;
332  sol = ope.inverse(rhs) ;
333  for (int i=0; i<nr; i++)
334  sol_hom1.set(lz, k, j, i) = sol(i) ;
335 
336  rhs.set(nr-2) = 0 ;
337  rhs.set(nr-1) = 1 ;
338  sol = ope.inverse(rhs) ;
339  for (int i=0; i<nr; i++)
340  sol_hom2.set(lz, k, j, i) = sol(i) ;
341 
342  }
343  }
344 
345  }
346  { int lz = nz-1 ;
347  double alpha = (*map).get_alpha()[lz] ;
348  for (int k=0; k < np; k++)
349  for (int j=0; j<nt; j++) {
350  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
351  if (nullite_plm(j, nt, k, np, base) == 1) {
352 
353  Matrice ope(nr,nr) ;
354  ope.annule_hard() ;
355  Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
356  Diff_sx2 sx2(base_r, nr) ; const Matrice& ms2 = sx2.get_matrice() ;
357 
358  ope = (mdx2 - l_q*(l_q+1)*ms2 + 2*l_q*ms2)/(alpha*alpha) ;
359 
360  for (int i=0; i<nr; i++)
361  ope.set(nr-1, i) = 0 ;
362  ope.set(nr-1, 0) = 1 ; //for the true homogeneous solution
363 
364  for (int i=0; i<nr; i++) {
365  ope.set(nr-2, i) = 1 ; //for the limit at inifinity
366  }
367 
368  if (l_q > 0) {
369  for (int i=0; i<nr; i++) {
370  ope.set(nr-3, i) = i*i ; //for the finite part (derivative = 0 at infty)
371  }
372  }
373  // cout << "l: " << l_q << endl ;
374 
375  Tbl rhs(nr) ;
376  rhs.annule_hard() ;
377  for (int i=0; i<nr; i++)
378  rhs.set(i) = (*source.get_spectral_va().c_cf)(lz, k, j, i) ;
379  if (l_q>0) rhs.set(nr-3) = 0 ;
380  rhs.set(nr-2) = 0 ;
381  rhs.set(nr-1) = 0 ;
382  ope.set_lu() ;
383  Tbl sol = ope.inverse(rhs) ;
384 
385  for (int i=0; i<nr; i++)
386  sol_part.set(lz, k, j, i) = sol(i) ;
387 
388  rhs.annule_hard() ;
389  rhs.set(nr-1) = 1 ;
390  sol = ope.inverse(rhs) ;
391  for (int i=0; i<nr; i++)
392  sol_hom2.set(lz, k, j, i) = sol(i) ;
393 
394  }
395  }
396  }
397 
398  Mtbl_cf dpart = sol_part ; dpart.dsdx() ;
399  Mtbl_cf dhom1 = sol_hom1 ; dhom1.dsdx() ;
400  Mtbl_cf dhom2 = sol_hom2 ; dhom2.dsdx() ;
401 
402 
403  // Now matching the homogeneous solutions between the different domains...
404 
405  for (int k=0; k < np; k++)
406  for (int j=0; j<nt; j++) {
407  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
408  if (nullite_plm(j, nt, k, np, base) == 1) {
409  Matrice systeme(2*nz-4, 2*nz-4) ;
410  systeme.annule_hard() ;
411  Tbl rhs(2*nz-4) ;
412  rhs.annule_hard() ;
413 
414  //First shell
415  int lin = 0 ;
416  int col = 0 ;
417 
418  double alpha = (*map).get_alpha()[1] ;
419 
420  systeme.set(lin, col) += sol_hom1.val_out_bound_jk(1, j, k) ;
421  rhs.set(lin) -= sol_part.val_out_bound_jk(1, j, k) ;
422 
423  lin++ ;
424  systeme.set(lin, col) += dhom1.val_out_bound_jk(1, j, k) / alpha ;
425  rhs.set(lin) -= dpart.val_out_bound_jk(1, j, k) / alpha ;
426  col += 1 ;
427 
428  //Shells
429  for (int lz=2; lz<nz-1; lz++) {
430  alpha = (*map).get_alpha()[lz] ;
431  lin-- ;
432  systeme.set(lin,col) -= sol_hom1.val_in_bound_jk(lz, j, k) ;
433  systeme.set(lin,col+1) -= sol_hom2.val_in_bound_jk(lz, j, k) ;
434  rhs.set(lin) += sol_part.val_in_bound_jk(lz, j, k) ;
435 
436  lin++ ;
437  systeme.set(lin,col) -= dhom1.val_in_bound_jk(lz, j, k) / alpha ;
438  systeme.set(lin,col+1) -= dhom2.val_in_bound_jk(lz, j, k) / alpha ;
439  rhs.set(lin) += dpart.val_in_bound_jk(lz, j, k) / alpha;
440 
441  lin++ ;
442  systeme.set(lin, col) += sol_hom1.val_out_bound_jk(lz, j, k) ;
443  systeme.set(lin, col+1) += sol_hom2.val_out_bound_jk(lz, j, k) ;
444  rhs.set(lin) -= sol_part.val_out_bound_jk(lz, j, k) ;
445 
446  lin++ ;
447  systeme.set(lin, col) += dhom1.val_out_bound_jk(lz, j, k) / alpha ;
448  systeme.set(lin, col+1) += dhom2.val_out_bound_jk(lz, j, k) / alpha ;
449  rhs.set(lin) -= dpart.val_out_bound_jk(lz, j, k) / alpha ;
450  col += 2 ;
451  }
452 
453  //CED
454  alpha = (*map).get_alpha()[nz-1] ;
455  lin-- ;
456  systeme.set(lin,col) -= sol_hom2.val_in_bound_jk(nz-1, j, k) ;
457  rhs.set(lin) += sol_part.val_in_bound_jk(nz-1, j, k) ;
458 
459  lin++ ;
460  systeme.set(lin,col) -= (-4*alpha)*dhom2.val_in_bound_jk(nz-1, j, k) ;
461  rhs.set(lin) += (-4*alpha)*dpart.val_in_bound_jk(nz-1, j, k) ;
462 
463  systeme.set_lu() ;
464 
465  // cout << systeme << endl;
466 
467  Tbl coef = systeme.inverse(rhs);
468  int indice = 0 ;
469 
470  // int tryr; cin >> tryr;
471  // for (int i=0; i<mgrid.get_nr(0); i++)
472  // sol_coef.set(0, k, j, i) = 0 ;
473  // sol_coef.set(0, k, j, 0) = (*bound.get_spectral_va().c_cf)(0, k, j, 0) ;
474 
475  for (int i=0; i<(*mgrid).get_nr(1); i++)
476  sol_coef.set(1, k, j, i) = sol_part(1, k, j, i)
477  +coef(indice)*sol_hom1(1, k, j, i) ;
478  indice +=1;
479 
480 
481  for (int lz=2; lz<nz-1; lz++) {
482  for (int i=0; i<(*mgrid).get_nr(lz); i++)
483  sol_coef.set(lz, k, j, i) = sol_part(lz, k, j, i)
484  +coef(indice)*sol_hom1(lz, k, j, i)
485  +coef(indice+1)*sol_hom2(lz, k, j, i) ;
486  indice += 2 ;
487  }
488  for (int i=0; i<(*mgrid).get_nr(nz-1); i++)
489  sol_coef.set(nz-1, k, j, i) = sol_part(nz-1, k, j, i)
490  +coef(indice)*sol_hom2(nz-1, k, j, i) ;
491 
492  }
493  }
494 
495 
496  delete phi.set_spectral_va().c ;
497  phi.set_spectral_va().c = 0x0 ;
498  // phi.set_spectral_va().ylm_i() ;
499 
500  phi.annule_domain(nz-1);
501 
502  resu = phi;
503 
504 }
505 }
Lorene
Lorene prototypes.
Definition: app_hor.h:67