5#include "utilitaires.h"
12 void tensorelliptic(
Scalar source,
Scalar& resu,
double fit,
double fit2,
double fit0d2,
double fit1d2,
double fit0d3,
double fit1d3) {
14 const int nz = (*source.get_mp().
get_mg()).get_nzone();
15 int nr = (*source.get_mp().
get_mg()).get_nr(1);
16 int nt = (*source.get_mp().
get_mg()).get_nt(1);
17 int np = (*source.get_mp().
get_mg()).get_np(1);
18 const Map_af* map =
dynamic_cast<const Map_af*
>(&source.get_mp()) ;
19 const Mg3d* mgrid = (*map).get_mg();
22 const Coord& rr = (*map).r;
25 rrr.set_spectral_va().set_base(source.get_spectral_va().base);
27 Scalar source_coq = source ;
30 source.set_spectral_va().ylm() ;
31 source_coq.set_spectral_va().ylm() ;
35 phi.set_spectral_va().set_base(source.get_spectral_va().base) ;
36 phi.set_spectral_va().ylm() ;
37 Mtbl_cf& sol_coef = (*
phi.set_spectral_va().c_cf) ;
39 const Base_val& base = source.get_spectral_base() ;
40 Mtbl_cf sol_part(mgrid, base) ; sol_part.annule_hard() ;
41 Mtbl_cf sol_hom1(mgrid, base) ; sol_hom1.annule_hard() ;
42 Mtbl_cf sol_hom2(mgrid, base) ; sol_hom2.annule_hard() ;
44 int l_q, m_q, base_r ;
47 for (
int k=0; k < np; k++)
48 for (
int j=0; j<nt; j++) {
49 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
50 if (nullite_plm(j, nt, k, np, base) == 1) {
51 for (
int ii=0 ; ii<nr ; ii++){
52 sol_hom1.set(lz, k, j, ii) = 0 ;
53 sol_part.set(lz, k, j, ii) = 0 ;
62 double alpha = (*map).get_alpha()[lz] ;
63 double beta = (*map).get_beta()[lz] ;
64 double ech = beta / alpha ;
65 for (
int k=0; k < np; k++)
66 for (
int j=0; j<nt; j++) {
67 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
68 if (nullite_plm(j, nt, k, np, base) == 1) {
85 ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
86 - 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)
87 + (fit)*alpha*(mx3dx2 + 2*ech*mx2dx2 + ech*ech*mxdx2) + fit2*alpha*alpha*(mx4dx2 + 2*ech*mx3dx2 + ech*ech*mx2dx2)
88 + 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));
93 for (
int col=0; col<nr; col++)
94 ope.set(nr-1, col) = 0 ;
95 ope.set(nr-1, 0) = 1 ;
100 for (
int i=0; i<nr; i++)
101 rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(1, k, j, i) ;
104 Tbl sol = ope.inverse(rhs) ;
107 for (
int i=0; i<nr; i++)
108 sol_part.
set(1, k, j, i) = sol(i) ;
112 sol = ope.inverse(rhs) ;
115 for (
int i=0; i<nr; i++)
116 sol_hom1.set(1, k, j, i) = sol(i) ;
129 double alpha = (*map).get_alpha()[lz] ;
130 double beta = (*map).get_beta()[lz] ;
131 double ech = beta / alpha ;
133 for (
int k=0; k < np; k++)
134 for (
int j=0; j<nt; j++) {
135 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
136 if (nullite_plm(j, nt, k, np, base) == 1) {
149 ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
150 - l_q*(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) ;
154 for (
int col=0; col<nr; col++)
155 ope.
set(nr-1, col) = 0 ;
156 ope.set(nr-1, 0) = 1 ;
157 for (
int col=0; col<nr; col++) {
158 ope.set(nr-2, col) = 0 ;
160 ope.set(nr-2, 1) = 1 ;
164 for (
int i=0; i<nr; i++)
165 rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(lz, k, j, i) ;
169 Tbl sol = ope.inverse(rhs) ;
171 for (
int i=0; i<nr; i++)
172 sol_part.
set(lz, k, j, i) = sol(i) ;
176 sol = ope.inverse(rhs) ;
177 for (
int i=0; i<nr; i++)
178 sol_hom1.set(lz, k, j, i) = sol(i) ;
182 sol = ope.inverse(rhs) ;
183 for (
int i=0; i<nr; i++)
184 sol_hom2.set(lz, k, j, i) = sol(i) ;
195 double alpha = (*map).get_alpha()[lz] ;
196 double beta = (*map).get_beta()[lz] ;
197 double ech = beta / alpha ;
199 for (
int k=0; k < np; k++)
200 for (
int j=0; j<nt; j++) {
201 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
202 if (nullite_plm(j, nt, k, np, base) == 1) {
215 ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
216 - l_q*(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) ;
218 for (
int col=0; col<nr; col++)
219 ope.
set(nr-1, col) = 0 ;
220 ope.set(nr-1, 0) = 1 ;
221 for (
int col=0; col<nr; col++) {
222 ope.set(nr-2, col) = 0 ;
224 ope.set(nr-2, 1) = 1 ;
228 for (
int i=0; i<nr; i++)
229 rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(lz, k, j, i) ;
233 Tbl sol = ope.inverse(rhs) ;
235 for (
int i=0; i<nr; i++)
236 sol_part.
set(lz, k, j, i) = sol(i) ;
240 sol = ope.inverse(rhs) ;
241 for (
int i=0; i<nr; i++)
242 sol_hom1.set(lz, k, j, i) = sol(i) ;
246 sol = ope.inverse(rhs) ;
247 for (
int i=0; i<nr; i++)
248 sol_hom2.set(lz, k, j, i) = sol(i) ;
258 for (
int lz=4; lz<nz-1; lz++) {
259 double ech = (*map).get_beta()[lz] / (*map).get_alpha()[lz] ;
260 for (
int k=0; k < np; k++)
261 for (
int j=0; j<nt; j++) {
262 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
263 if (nullite_plm(j, nt, k, np, base) == 1) {
274 ope = mx2dx2 + 2*ech*mxdx2 + ech*ech*mdx2 + 2*(mxdx + ech*mdx)
277 for (
int col=0; col<nr; col++)
278 ope.
set(nr-1, col) = 0 ;
279 ope.set(nr-1, 0) = 1 ;
280 for (
int col=0; col<nr; col++) {
281 ope.set(nr-2, col) = 0 ;
283 ope.set(nr-2, 1) = 1 ;
287 for (
int i=0; i<nr; i++)
288 rhs.set(i) = (*source_coq.get_spectral_va().c_cf)(lz, k, j, i) ;
292 Tbl sol = ope.inverse(rhs) ;
294 for (
int i=0; i<nr; i++)
295 sol_part.
set(lz, k, j, i) = sol(i) ;
299 sol = ope.inverse(rhs) ;
300 for (
int i=0; i<nr; i++)
301 sol_hom1.set(lz, k, j, i) = sol(i) ;
305 sol = ope.inverse(rhs) ;
306 for (
int i=0; i<nr; i++)
307 sol_hom2.set(lz, k, j, i) = sol(i) ;
314 double alpha = (*map).get_alpha()[lz] ;
315 for (
int k=0; k < np; k++)
316 for (
int j=0; j<nt; j++) {
317 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
318 if (nullite_plm(j, nt, k, np, base) == 1) {
325 ope = (mdx2 - l_q*(l_q+1)*ms2)/(alpha*alpha) ;
327 for (
int i=0; i<nr; i++)
328 ope.
set(nr-1, i) = 0 ;
329 ope.set(nr-1, 0) = 1 ;
331 for (
int i=0; i<nr; i++) {
332 ope.set(nr-2, i) = 1 ;
336 for (
int i=0; i<nr; i++) {
337 ope.set(nr-3, i) = i*i ;
343 for (
int i=0; i<nr; i++)
344 rhs.set(i) = (*source.get_spectral_va().c_cf)(lz, k, j, i) ;
345 if (l_q>0) rhs.set(nr-3) = 0 ;
349 Tbl sol = ope.inverse(rhs) ;
351 for (
int i=0; i<nr; i++)
352 sol_part.
set(lz, k, j, i) = sol(i) ;
356 sol = ope.inverse(rhs) ;
357 for (
int i=0; i<nr; i++)
358 sol_hom2.set(lz, k, j, i) = sol(i) ;
371 for (
int k=0; k < np; k++)
372 for (
int j=0; j<nt; j++) {
373 base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
374 if (nullite_plm(j, nt, k, np, base) == 1) {
375 Matrice systeme(2*nz-4, 2*nz-4) ;
376 systeme.annule_hard() ;
384 double alpha = (*map).get_alpha()[1] ;
386 systeme.set(lin, col) += sol_hom1.val_out_bound_jk(1, j, k) ;
387 rhs.set(lin) -= sol_part.val_out_bound_jk(1, j, k) ;
390 systeme.set(lin, col) += dhom1.val_out_bound_jk(1, j, k) / alpha ;
391 rhs.set(lin) -= dpart.val_out_bound_jk(1, j, k) / alpha ;
395 for (
int lz=2; lz<nz-1; lz++) {
396 alpha = (*map).get_alpha()[lz] ;
398 systeme.set(lin,col) -= sol_hom1.val_in_bound_jk(lz, j, k) ;
399 systeme.set(lin,col+1) -= sol_hom2.val_in_bound_jk(lz, j, k) ;
400 rhs.set(lin) += sol_part.val_in_bound_jk(lz, j, k) ;
403 systeme.set(lin,col) -= dhom1.val_in_bound_jk(lz, j, k) / alpha ;
404 systeme.set(lin,col+1) -= dhom2.val_in_bound_jk(lz, j, k) / alpha ;
405 rhs.set(lin) += dpart.val_in_bound_jk(lz, j, k) / alpha;
408 systeme.set(lin, col) += sol_hom1.val_out_bound_jk(lz, j, k) ;
409 systeme.set(lin, col+1) += sol_hom2.val_out_bound_jk(lz, j, k) ;
410 rhs.set(lin) -= sol_part.val_out_bound_jk(lz, j, k) ;
413 systeme.set(lin, col) += dhom1.val_out_bound_jk(lz, j, k) / alpha ;
414 systeme.set(lin, col+1) += dhom2.val_out_bound_jk(lz, j, k) / alpha ;
415 rhs.set(lin) -= dpart.val_out_bound_jk(lz, j, k) / alpha ;
420 alpha = (*map).get_alpha()[nz-1] ;
422 systeme.set(lin,col) -= sol_hom2.val_in_bound_jk(nz-1, j, k) ;
423 rhs.set(lin) += sol_part.val_in_bound_jk(nz-1, j, k) ;
426 systeme.set(lin,col) -= (-4*alpha)*dhom2.val_in_bound_jk(nz-1, j, k) ;
427 rhs.set(lin) += (-4*alpha)*dpart.val_in_bound_jk(nz-1, j, k) ;
431 Tbl coef = systeme.inverse(rhs);
435 for (
int i=0; i<(*mgrid).get_nr(1); i++)
436 sol_coef.
set(1, k, j, i) = sol_part(1, k, j, i)
437 +coef(indice)*sol_hom1(1, k, j, i) ;
441 for (
int lz=2; lz<nz-1; lz++) {
442 for (
int i=0; i<(*mgrid).get_nr(lz); i++)
443 sol_coef.set(lz, k, j, i) = sol_part(lz, k, j, i)
444 +coef(indice)*sol_hom1(lz, k, j, i)
445 +coef(indice+1)*sol_hom2(lz, k, j, i) ;
448 for (
int i=0; i<(*mgrid).get_nr(nz-1); i++)
449 sol_coef.set(nz-1, k, j, i) = sol_part(nz-1, k, j, i)
450 +coef(indice)*sol_hom2(nz-1, k, j, i) ;
456 delete phi.set_spectral_va().c ;
457 phi.set_spectral_va().c = 0x0 ;
458 phi.set_spectral_va().ylm_i() ;
460 phi.annule_domain(nz-1);
Bases of the spectral expansions.
Active physical coordinates and mapping derivatives.
Class for the elementary differential operator (see the base class Diff ).
Class for the elementary differential operator (see the base class Diff ).
Class for the elementary differential operator Identity (see the base class Diff ).
Class for the elementary differential operator division by (see the base class Diff ).
Class for the elementary differential operator (see the base class Diff ).
Class for the elementary differential operator (see the base class Diff ).
Class for the elementary differential operator (see the base class Diff ).
Class for the elementary differential operator (see the base class Diff ).
Class for the elementary differential operator (see the base class Diff ).
double & set(int j, int i)
Read/write of a particuliar element.
Coefficients storage for the multi-domain spectral method.
Tensor field of valence 0 (or component of a tensorial field).
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
double & set(int i)
Read/write of a particular element (index i) (1D case).
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Coord phi
coordinate centered on the grid