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
11namespace 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::Base_val
Bases of the spectral expansions.
Definition base_val.h:325
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:2042
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:279
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:196
Lorene::Mtbl_cf::dsdx
void dsdx()
Definition valeur_dsdx.C:156
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:393
Lorene::Scalar::mult_r
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
Definition scalar_r_manip.C:211
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:67
Lorene::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:777
Lorene::phi
Coord phi
coordinate centered on the grid
Definition map.h:732