LORENE
tensorellipticCt.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 , associated with spectral variable tilde (C) with main characteristics
13 //fit and fit2 (Suitable for tensorial resolution of BH spacetime)
14 // See also tilde(laplacian)
15 
16  void tensorellipticCt( 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 
28  // Some helpful stuff...
29 
30  const Coord& rr = (*map).r ;
31  Scalar rrr (*map) ;
32  rrr = rr ;
33  rrr.set_spectral_va().set_base(source.get_spectral_va().base);
34 
35  Scalar source_coq = source ;
36  source_coq.mult_r() ;
37  source_coq.mult_r() ;
38  source.set_spectral_va().ylm() ;
39  source_coq.set_spectral_va().ylm() ;
40  Scalar phi(source.get_mp()) ;
41  phi.annule_hard() ;
42  // phi.std_spectral_base();
43  phi.set_spectral_va().set_base(source.get_spectral_va().base) ;
44  phi.set_spectral_va().ylm() ;
45  Mtbl_cf& sol_coef = (*phi.set_spectral_va().c_cf) ;
46 
47  const Base_val& base = source.get_spectral_base() ;
48  Mtbl_cf sol_part(mgrid, base) ; sol_part.annule_hard() ;
49  Mtbl_cf sol_hom1(mgrid, base) ; sol_hom1.annule_hard() ;
50  Mtbl_cf sol_hom2(mgrid, base) ; sol_hom2.annule_hard() ;
51 
52  int l_q, m_q, base_r ;
53 
54  { int lz = 0 ;
55  for (int k=0; k < np; k++)
56  for (int j=0; j<nt; j++) {
57  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
58  if (nullite_plm(j, nt, k, np, base) == 1) {
59  for (int ii=0 ; ii<nr ; ii++){
60  sol_hom1.set(lz, k, j, ii) = 0 ;
61  sol_part.set(lz, k, j, ii) = 0 ;
62  }
63 
64  }
65  }
66  }
67 
68 
69  { 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.
70  double alpha = (*map).get_alpha()[lz] ;
71  double beta = (*map).get_beta()[lz] ;
72  double ech = beta / alpha ;
73  for (int k=0; k < np; k++)
74  for (int j=0; j<nt; j++) {
75  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
76  if (nullite_plm(j, nt, k, np, base) == 1) {
77 
78 
79 
80  Matrice ope(nr,nr) ;
81  ope.annule_hard() ;
82 
83  Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
84  Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
85  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
86  Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
87  Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
88  Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
89  Diff_x3dsdx2 x3dx2 (base_r, nr); const Matrice& mx3dx2 = x3dx2.get_matrice();
90  Diff_x4dsdx2 x4dx2 (base_r, nr); const Matrice& mx4dx2 = x4dx2.get_matrice();
91 
92 
93 
94 // ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
95 // - 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) ;
96 
97 
98 // ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
99 // - 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)
100 // + (fit)*alpha*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2) + fit2*alpha*alpha*(mx4dx2 + 2*ech*mx3dx2 + ech*ech*mx2dx2)
101 // + 2*fit2*alpha*(beta-1.)*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2) + fit2*(beta-1.)*(beta-1.)*(mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2));
102 
103 
104 
105  ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
106  - l_q*(l_q+1)*mid -2*(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)
107  + (fit)*alpha*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2) + fit2*alpha*alpha*(mx4dx2 + 2*ech*mx3dx2 + ech*ech*mx2dx2)
108  + 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));
109 
110 
111 
112 
113  for (int col=0; col<nr; col++)
114  ope.set(nr-1, col) = 0 ;
115  ope.set(nr-1, 0) = 1 ;
116 
117 
118  Tbl rhs(nr);
119  rhs.annule_hard() ;
120  for (int i=0; i<nr; i++)
121  rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(1, k, j, i) ;
122  rhs.set(nr-1) = 0 ;
123  ope.set_lu() ;
124  Tbl sol = ope.inverse(rhs) ;
125 
126 
127  for (int i=0; i<nr; i++)
128  sol_part.set(1, k, j, i) = sol(i) ;
129 
130  rhs.annule_hard();
131  rhs.set(nr-1) = 1 ;
132  sol = ope.inverse(rhs) ;
133 
134 
135  for (int i=0; i<nr; i++)
136  sol_hom1.set(1, k, j, i) = sol(i) ;
137 
138 
139  }
140  }
141  }
142 
143  // Attention! zones 2 et 3traitee separement egalement!!!
144 
145 
146  // Current implementations only allow grids with more than 3 shells.
147 
148  { int lz = 2 ;
149 
150  double alpha = (*map).get_alpha()[lz] ;
151  double beta = (*map).get_beta()[lz] ;
152  double ech = beta / alpha ;
153 
154  for (int k=0; k < np; k++)
155  for (int j=0; j<nt; j++) {
156  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
157  if (nullite_plm(j, nt, k, np, base) == 1) {
158 
159  Matrice ope(nr,nr) ;
160  ope.annule_hard() ;
161 
162  Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
163  Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
164  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
165  Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
166  Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
167  Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
168  Diff_x3dsdx2 x3dx2 (base_r, nr); const Matrice& mx3dx2 = x3dx2.get_matrice();
169 
170  ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
171  - l_q*(l_q+1)*mid -2*(l_q+1)*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) ;
172 
173 
174 
175 // ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
176 // - l_q*(l_q+1)*mid ;
177 
178  for (int col=0; col<nr; col++)
179  ope.set(nr-1, col) = 0 ;
180  ope.set(nr-1, 0) = 1 ;
181  for (int col=0; col<nr; col++) {
182  ope.set(nr-2, col) = 0 ;
183  }
184  ope.set(nr-2, 1) = 1 ;
185 
186  Tbl rhs(nr) ;
187  rhs.annule_hard() ;
188  for (int i=0; i<nr; i++)
189  rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(lz, k, j, i) ;
190  rhs.set(nr-2) = 0 ;
191  rhs.set(nr-1) = 0 ;
192  ope.set_lu() ;
193  Tbl sol = ope.inverse(rhs) ;
194 
195  for (int i=0; i<nr; i++)
196  sol_part.set(lz, k, j, i) = sol(i) ;
197 
198  rhs.annule_hard() ;
199  rhs.set(nr-2) = 1 ;
200  sol = ope.inverse(rhs) ;
201  for (int i=0; i<nr; i++)
202  sol_hom1.set(lz, k, j, i) = sol(i) ;
203 
204  rhs.set(nr-2) = 0 ;
205  rhs.set(nr-1) = 1 ;
206  sol = ope.inverse(rhs) ;
207  for (int i=0; i<nr; i++)
208  sol_hom2.set(lz, k, j, i) = sol(i) ;
209 
210  }
211  }
212 
213  }
214 
215 
216 
217  { int lz = 3 ;
218 
219  double alpha = (*map).get_alpha()[lz] ;
220  double beta = (*map).get_beta()[lz] ;
221  double ech = beta / alpha ;
222 
223  for (int k=0; k < np; k++)
224  for (int j=0; j<nt; j++) {
225  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
226  if (nullite_plm(j, nt, k, np, base) == 1) {
227 
228  Matrice ope(nr,nr) ;
229  ope.annule_hard() ;
230 
231  Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
232  Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
233  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
234  Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
235  Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
236  Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
237  Diff_x3dsdx2 x3dx2 (base_r, nr); const Matrice& mx3dx2 = x3dx2.get_matrice();
238 
239  ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
240  - l_q*(l_q+1)*mid -2*(l_q+1)*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) ;
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  // Current implementations only allow grids with more than 2 shells.
286 
287  for (int lz=4; lz<nz-1; lz++) {
288  double ech = (*map).get_beta()[lz] / (*map).get_alpha()[lz] ;
289  for (int k=0; k < np; k++)
290  for (int j=0; j<nt; j++) {
291  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
292  if (nullite_plm(j, nt, k, np, base) == 1) {
293 
294  Matrice ope(nr,nr) ;
295  ope.annule_hard() ;
296 
297  Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
298  Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
299  Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
300  Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
301  Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
302  Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
303  ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
304  - l_q*(l_q+1)*mid -2*(l_q+1)*mid;
305 
306  for (int col=0; col<nr; col++)
307  ope.set(nr-1, col) = 0 ;
308  ope.set(nr-1, 0) = 1 ;
309  for (int col=0; col<nr; col++) {
310  ope.set(nr-2, col) = 0 ;
311  }
312  ope.set(nr-2, 1) = 1 ;
313 
314  Tbl rhs(nr) ;
315  rhs.annule_hard() ;
316  for (int i=0; i<nr; i++)
317  rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(lz, k, j, i) ;
318  rhs.set(nr-2) = 0 ;
319  rhs.set(nr-1) = 0 ;
320  ope.set_lu() ;
321  Tbl sol = ope.inverse(rhs) ;
322 
323  for (int i=0; i<nr; i++)
324  sol_part.set(lz, k, j, i) = sol(i) ;
325 
326  rhs.annule_hard() ;
327  rhs.set(nr-2) = 1 ;
328  sol = ope.inverse(rhs) ;
329  for (int i=0; i<nr; i++)
330  sol_hom1.set(lz, k, j, i) = sol(i) ;
331 
332  rhs.set(nr-2) = 0 ;
333  rhs.set(nr-1) = 1 ;
334  sol = ope.inverse(rhs) ;
335  for (int i=0; i<nr; i++)
336  sol_hom2.set(lz, k, j, i) = sol(i) ;
337 
338  }
339  }
340 
341  }
342  { int lz = nz-1 ;
343  double alpha = (*map).get_alpha()[lz] ;
344  for (int k=0; k < np; k++)
345  for (int j=0; j<nt; j++) {
346  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
347  if (nullite_plm(j, nt, k, np, base) == 1) {
348 
349  Matrice ope(nr,nr) ;
350  ope.annule_hard() ;
351  Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
352  Diff_sx2 sx2(base_r, nr) ; const Matrice& ms2 = sx2.get_matrice() ;
353 
354  ope = (mdx2 - l_q*(l_q+1)*ms2 -2*(l_q+1)*ms2)/(alpha*alpha) ;
355 
356  for (int i=0; i<nr; i++)
357  ope.set(nr-1, i) = 0 ;
358  ope.set(nr-1, 0) = 1 ; //for the true homogeneous solution
359 
360  for (int i=0; i<nr; i++) {
361  ope.set(nr-2, i) = 1 ; //for the limit at inifinity
362  }
363 
364  if (l_q > 0) {
365  for (int i=0; i<nr; i++) {
366  ope.set(nr-3, i) = i*i ; //for the finite part (derivative = 0 at infty)
367  }
368  }
369  // cout << "l: " << l_q << endl ;
370 
371  Tbl rhs(nr) ;
372  rhs.annule_hard() ;
373  for (int i=0; i<nr; i++)
374  rhs.set(i) = (*source.get_spectral_va().c_cf)(lz, k, j, i) ;
375  if (l_q>0) rhs.set(nr-3) = 0 ;
376  rhs.set(nr-2) = 0 ;
377  rhs.set(nr-1) = 0 ;
378  ope.set_lu() ;
379  Tbl sol = ope.inverse(rhs) ;
380 
381  for (int i=0; i<nr; i++)
382  sol_part.set(lz, k, j, i) = sol(i) ;
383 
384  rhs.annule_hard() ;
385  rhs.set(nr-1) = 1 ;
386  sol = ope.inverse(rhs) ;
387  for (int i=0; i<nr; i++)
388  sol_hom2.set(lz, k, j, i) = sol(i) ;
389 
390  }
391  }
392  }
393 
394  Mtbl_cf dpart = sol_part ; dpart.dsdx() ;
395  Mtbl_cf dhom1 = sol_hom1 ; dhom1.dsdx() ;
396  Mtbl_cf dhom2 = sol_hom2 ; dhom2.dsdx() ;
397 
398 
399 // cout << sol_hom1 << endl;
400 // int thj; cin >> thj;
401 // cout << sol_hom2 << endl;
402 // int thj2; cin >> thj2;
403 // cout << sol_part << endl;
404 // int thj3; cin >> thj3;
405 
406 
407 
408 
409  // Now matching the homogeneous solutions between the different domains...
410 
411  for (int k=0; k < np; k++)
412  for (int j=0; j<nt; j++) {
413  base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
414  if (nullite_plm(j, nt, k, np, base) == 1) {
415  Matrice systeme(2*nz-4, 2*nz-4) ;
416  systeme.annule_hard() ;
417  Tbl rhs(2*nz-4) ;
418  rhs.annule_hard() ;
419 
420  //First shell
421  int lin = 0 ;
422  int col = 0 ;
423 
424  double alpha = (*map).get_alpha()[1] ;
425 
426  systeme.set(lin, col) += sol_hom1.val_out_bound_jk(1, j, k) ;
427  rhs.set(lin) -= sol_part.val_out_bound_jk(1, j, k) ;
428 
429  lin++ ;
430  systeme.set(lin, col) += dhom1.val_out_bound_jk(1, j, k) / alpha ;
431  rhs.set(lin) -= dpart.val_out_bound_jk(1, j, k) / alpha ;
432  col += 1 ;
433 
434  //Shells
435  for (int lz=2; lz<nz-1; lz++) {
436  alpha = (*map).get_alpha()[lz] ;
437  lin-- ;
438  systeme.set(lin,col) -= sol_hom1.val_in_bound_jk(lz, j, k) ;
439  systeme.set(lin,col+1) -= sol_hom2.val_in_bound_jk(lz, j, k) ;
440  rhs.set(lin) += sol_part.val_in_bound_jk(lz, j, k) ;
441 
442  lin++ ;
443  systeme.set(lin,col) -= dhom1.val_in_bound_jk(lz, j, k) / alpha ;
444  systeme.set(lin,col+1) -= dhom2.val_in_bound_jk(lz, j, k) / alpha ;
445  rhs.set(lin) += dpart.val_in_bound_jk(lz, j, k) / alpha;
446 
447  lin++ ;
448  systeme.set(lin, col) += sol_hom1.val_out_bound_jk(lz, j, k) ;
449  systeme.set(lin, col+1) += sol_hom2.val_out_bound_jk(lz, j, k) ;
450  rhs.set(lin) -= sol_part.val_out_bound_jk(lz, j, k) ;
451 
452  lin++ ;
453  systeme.set(lin, col) += dhom1.val_out_bound_jk(lz, j, k) / alpha ;
454  systeme.set(lin, col+1) += dhom2.val_out_bound_jk(lz, j, k) / alpha ;
455  rhs.set(lin) -= dpart.val_out_bound_jk(lz, j, k) / alpha ;
456  col += 2 ;
457  }
458 
459  //CED
460  alpha = (*map).get_alpha()[nz-1] ;
461  lin-- ;
462  systeme.set(lin,col) -= sol_hom2.val_in_bound_jk(nz-1, j, k) ;
463  rhs.set(lin) += sol_part.val_in_bound_jk(nz-1, j, k) ;
464 
465  lin++ ;
466  systeme.set(lin,col) -= (-4*alpha)*dhom2.val_in_bound_jk(nz-1, j, k) ;
467  rhs.set(lin) += (-4*alpha)*dpart.val_in_bound_jk(nz-1, j, k) ;
468 
469  systeme.set_lu() ;
470 
471  // cout << systeme << endl;
472 
473  Tbl coef = systeme.inverse(rhs);
474  int indice = 0 ;
475 
476  // int tryr; cin >> tryr;
477  // for (int i=0; i<mgrid.get_nr(0); i++)
478  // sol_coef.set(0, k, j, i) = 0 ;
479  // sol_coef.set(0, k, j, 0) = (*bound.get_spectral_va().c_cf)(0, k, j, 0) ;
480 
481  for (int i=0; i<(*mgrid).get_nr(1); i++)
482  sol_coef.set(1, k, j, i) = sol_part(1, k, j, i)
483  +coef(indice)*sol_hom1(1, k, j, i) ;
484  indice +=1;
485 
486 
487  for (int lz=2; lz<nz-1; lz++) {
488  for (int i=0; i<(*mgrid).get_nr(lz); i++)
489  sol_coef.set(lz, k, j, i) = sol_part(lz, k, j, i)
490  +coef(indice)*sol_hom1(lz, k, j, i)
491  +coef(indice+1)*sol_hom2(lz, k, j, i) ;
492  indice += 2 ;
493  }
494  for (int i=0; i<(*mgrid).get_nr(nz-1); i++)
495  sol_coef.set(nz-1, k, j, i) = sol_part(nz-1, k, j, i)
496  +coef(indice)*sol_hom2(nz-1, k, j, i) ;
497 
498  }
499  }
500 
501 
502 
503  // cout << "sol coef" << endl;
504 
505 // int trg; cin >> trg;
506 // cout << sol_coef << endl;
507 
508 // int dn; cin >> dn;
509 
510  for (int lz=0; lz<nz; lz++) {
511  for (int i=0; i<(*mgrid).get_nr(lz); i++)
512 
513  sol_coef.set(lz,0,0,i) = 0.; // Probleme de définition de l'operateur pour k=j=0; a revoir.
514 
515  }
516 
517 
518 
519 // cout << "sol coef" << endl;
520 
521 // cin >> trg;
522 // cout << sol_coef << endl;
523 
524 // cin >> dn;
525 
526 
527 
528  delete phi.set_spectral_va().c ;
529  phi.set_spectral_va().c = 0x0 ;
530  // phi.set_spectral_va().ylm_i() ;
531  // phi.set_spectral_va().ylm();
532 
533  phi.annule_domain(nz-1);
534 
535  resu = phi;
536 
537 }
538 }
Lorene
Lorene prototypes.
Definition: app_hor.h:67