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
11namespace 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::Base_val
Bases of the spectral expansions.
Definition base_val.h:322
Lorene::Coord
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Lorene::Diff_dsdx2
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:172
Lorene::Diff_dsdx
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:129
Lorene::Diff_id
Class for the elementary differential operator Identity (see the base class Diff ).
Definition diff.h:210
Lorene::Diff_sx2
Class for the elementary differential operator division by (see the base class Diff ).
Definition diff.h:369
Lorene::Diff_x2dsdx2
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:490
Lorene::Diff_x3dsdx2
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:652
Lorene::Diff_x4dsdx2
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:692
Lorene::Diff_xdsdx2
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:531
Lorene::Diff_xdsdx
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:409
Lorene::Map_af
Affine radial mapping.
Definition map.h:2027
Lorene::Matrice
Matrix handling.
Definition matrice.h:152
Lorene::Matrice::set
double & set(int j, int i)
Read/write of a particuliar element.
Definition matrice.h:277
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Lorene::Mtbl_cf::dsdx
void dsdx()
Definition valeur_dsdx.C:153
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::mult_r
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
Definition scalar_r_manip.C:208
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Lorene::phi
Coord phi
coordinate centered on the grid
Definition map.h:720