26char pois_vect_r0_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/pois_vect_r0.C,v 1.2 2014/10/13 08:53:29 j_novak Exp $" ;
59Scalar pois_vect_r0(
const Scalar& source) {
61 const Map& map0 = source.
get_mp() ;
62 const Map_af* map1 =
dynamic_cast<const Map_af*
>(&map0) ;
64 const Map_af& map = *map1 ;
66 const Mg3d& mgrid = *map.get_mg() ;
67 int nz = mgrid.get_nzone() ;
70 if (source.get_etat() == ETATZERO) {
76 resu.std_spectral_base_odd() ;
77 resu.set_spectral_va().ylm() ;
78 Mtbl_cf& sol_coef = (*resu.set_spectral_va().c_cf) ;
79 const Base_val& base = source.get_spectral_base() ;
80 assert(resu.get_spectral_base() == base) ;
81 assert(source.check_dzpuis(4)) ;
83 Mtbl_cf sol_part(mgrid, base) ; sol_part.annule_hard() ;
84 Mtbl_cf sol_hom1(mgrid, base) ; sol_hom1.annule_hard() ;
85 Mtbl_cf sol_hom2(mgrid, base) ; sol_hom2.annule_hard() ;
88 int nr = mgrid.get_nr(lz) ;
89 double alpha2 = map.get_alpha()[lz]*map.get_alpha()[lz] ;
90 assert(mgrid.get_type_r(lz) == RARE) ;
94 Diff_dsdx2 dx2(base_r, nr) ;
const Matrice& mdx2 = dx2.get_matrice() ;
95 Diff_sxdsdx sdx(base_r, nr) ;
const Matrice& msdx = sdx.get_matrice() ;
96 Diff_sx2 sx2(base_r, nr) ;
const Matrice& ms2 = sx2.get_matrice() ;
98 for (
int lin=0; lin<nr-1; lin++)
99 for (
int col=0; col<nr; col++)
100 ope.set(lin,col) = (mdx2(lin,col) + 2*msdx(lin,col) - 2*ms2(lin,col))/alpha2 ;
102 ope.set(nr-1, 0) = 1 ;
103 for (
int i=1; i<nr; i++)
104 ope.set(nr-1, i) = 0 ;
108 for (
int i=0; i<nr-1; i++)
109 rhs.set(i) = (*source.get_spectral_va().c_cf)(lz, 0, 0, i) ;
112 Tbl sol = ope.inverse(rhs) ;
114 for (
int i=0; i<nr; i++)
115 sol_part.set(lz, 0, 0, i) = sol(i) ;
119 sol = ope.inverse(rhs) ;
120 for (
int i=0; i<nr; i++)
121 sol_hom1.set(lz, 0, 0, i) = sol(i) ;
125 for (
int lz=1; lz<nz-1; lz++) {
126 int nr = mgrid.get_nr(lz) ;
127 double alpha = map.get_alpha()[lz] ;
128 double beta = map.get_beta()[lz] ;
129 double ech = beta / alpha ;
130 assert(mgrid.get_type_r(lz) == FIN) ;
136 Diff_dsdx dx(base_r, nr) ;
const Matrice& mdx = dx.get_matrice() ;
137 Diff_dsdx2 dx2(base_r, nr) ;
const Matrice& mdx2 = dx2.get_matrice() ;
138 Diff_id id(base_r, nr) ;
const Matrice& mid =
id.get_matrice() ;
139 Diff_xdsdx xdx(base_r, nr) ;
const Matrice& mxdx = xdx.get_matrice() ;
140 Diff_xdsdx2 xdx2(base_r, nr) ;
const Matrice& mxdx2 = xdx2.get_matrice() ;
141 Diff_x2dsdx2 x2dx2(base_r, nr) ;
const Matrice& mx2dx2 = x2dx2.get_matrice() ;
143 for (
int lin=0; lin<nr-2; lin++)
144 for (
int col=0; col<nr; col++)
145 ope.set(lin, col) = mx2dx2(lin, col) + 2*ech*mxdx2(lin, col) + ech*ech*mdx2(lin, col)
146 + 2*(mxdx(lin, col) + ech*mdx(lin, col)) - 2*mid(lin, col) ;
148 ope.set(nr-2, 0) = 0 ;
149 ope.set(nr-2, 1) = 1 ;
150 for (
int col=2; col<nr; col++) {
151 ope.set(nr-2, col) = 0 ;
153 ope.set(nr-1, 0) = 1 ;
154 for (
int col=1; col<nr; col++)
155 ope.set(nr-1, col) = 0 ;
159 for (
int i=0; i<nr; i++)
160 src.set(i) = (*source.get_spectral_va().c_cf)(lz, 0, 0, i) ;
161 Tbl xsrc = src ; multx_1d(nr, &xsrc.t, base_r) ;
162 Tbl x2src = src ; multx2_1d(nr, &x2src.t, base_r) ;
165 for (
int i=0; i<nr-2; i++)
166 rhs.set(i) = beta*beta*src(i) + 2*beta*alpha*xsrc(i) + alpha*alpha*x2src(i) ;
170 Tbl sol = ope.inverse(rhs) ;
172 for (
int i=0; i<nr; i++)
173 sol_part.set(lz, 0, 0, i) = sol(i) ;
177 sol = ope.inverse(rhs) ;
178 for (
int i=0; i<nr; i++)
179 sol_hom1.set(lz, 0, 0, i) = sol(i) ;
183 sol = ope.inverse(rhs) ;
184 for (
int i=0; i<nr; i++)
185 sol_hom2.set(lz, 0, 0, i) = sol(i) ;
190 int nr = mgrid.get_nr(lz) ;
191 double alpha2 = map.get_alpha()[lz]*map.get_alpha()[lz] ;
192 assert(mgrid.get_type_r(lz) == UNSURR) ;
197 Diff_dsdx2 dx2(base_r, nr) ;
const Matrice& mdx2 = dx2.get_matrice() ;
198 Diff_sx2 sx2(base_r, nr) ;
const Matrice& ms2 = sx2.get_matrice() ;
200 for (
int lin=0; lin<nr-3; lin++)
201 for (
int col=0; col<nr; col++)
202 ope.set(lin, col) = (mdx2(lin, col) - 2*ms2(lin, col))/alpha2 ;
204 for (
int i=0; i<nr; i++) {
205 ope.set(nr-3, i) = i*i ;
208 for (
int i=0; i<nr; i++) {
209 ope.set(nr-2, i) = 1 ;
212 ope.set(nr-1, 0) = 1 ;
213 for (
int i=1; i<nr; i++)
214 ope.set(nr-1, i) = 0 ;
218 for (
int i=0; i<nr-3; i++)
219 rhs.set(i) = (*source.get_spectral_va().c_cf)(lz, 0, 0, i) ;
224 Tbl sol = ope.inverse(rhs) ;
225 for (
int i=0; i<nr; i++)
226 sol_part.set(lz, 0, 0, i) = sol(i) ;
230 sol = ope.inverse(rhs) ;
231 for (
int i=0; i<nr; i++)
232 sol_hom2.set(lz, 0, 0, i) = sol(i) ;
236 Mtbl_cf dpart = sol_part ; dpart.dsdx() ;
237 Mtbl_cf dhom1 = sol_hom1 ; dhom1.dsdx() ;
238 Mtbl_cf dhom2 = sol_hom2 ; dhom2.dsdx() ;
240 Matrice systeme(2*(nz-1), 2*(nz-1)) ;
241 systeme.annule_hard() ;
248 double alpha = map.get_alpha()[0] ;
249 systeme.set(lin,col) = sol_hom1.val_out_bound_jk(0, 0, 0) ;
250 rhs.set(lin) -= sol_part.val_out_bound_jk(0, 0, 0) ;
252 systeme.set(lin,col) = dhom1.val_out_bound_jk(0, 0, 0) / alpha ;
253 rhs.set(lin) -= dpart.val_out_bound_jk(0, 0, 0) / alpha ;
257 for (
int lz=1; lz<nz-1; lz++) {
258 alpha = map.get_alpha()[lz] ;
260 systeme.set(lin,col) -= sol_hom1.val_in_bound_jk(lz, 0, 0) ;
261 systeme.set(lin,col+1) -= sol_hom2.val_in_bound_jk(lz, 0, 0) ;
262 rhs.set(lin) += sol_part.val_in_bound_jk(lz, 0, 0) ;
265 systeme.set(lin,col) -= dhom1.val_in_bound_jk(lz, 0, 0) / alpha ;
266 systeme.set(lin,col+1) -= dhom2.val_in_bound_jk(lz, 0, 0) / alpha ;
267 rhs.set(lin) += dpart.val_in_bound_jk(lz, 0, 0) / alpha;
270 systeme.set(lin, col) += sol_hom1.val_out_bound_jk(lz, 0, 0) ;
271 systeme.set(lin, col+1) += sol_hom2.val_out_bound_jk(lz, 0, 0) ;
272 rhs.set(lin) -= sol_part.val_out_bound_jk(lz, 0, 0) ;
275 systeme.set(lin, col) += dhom1.val_out_bound_jk(lz, 0, 0) / alpha ;
276 systeme.set(lin, col+1) += dhom2.val_out_bound_jk(lz, 0, 0) / alpha ;
277 rhs.set(lin) -= dpart.val_out_bound_jk(lz, 0, 0) / alpha ;
282 alpha = map.get_alpha()[nz-1] ;
284 systeme.set(lin,col) -= sol_hom2.val_in_bound_jk(nz-1, 0, 0) ;
285 rhs.set(lin) += sol_part.val_in_bound_jk(nz-1, 0, 0) ;
288 systeme.set(lin,col) -= (-4*alpha)*dhom2.val_in_bound_jk(nz-1, 0, 0) ;
289 rhs.set(lin) += (-4*alpha)*dpart.val_in_bound_jk(nz-1, 0, 0) ;
292 Tbl coef = systeme.inverse(rhs) ;
295 for (
int i=0; i<mgrid.get_nr(0); i++)
296 sol_coef.set(0, 0, 0, i) = sol_part(0, 0, 0, i)
297 + coef(indice)*sol_hom1(0, 0, 0, i) ;
299 for (
int lz=1; lz<nz-1; lz++) {
300 for (
int i=0; i<mgrid.get_nr(lz); i++)
301 sol_coef.set(lz, 0, 0, i) = sol_part(lz, 0, 0, i)
302 +coef(indice)*sol_hom1(lz, 0, 0, i)
303 +coef(indice+1)*sol_hom2(lz, 0, 0, i) ;
306 for (
int i=0; i<mgrid.get_nr(nz-1); i++)
307 sol_coef.set(nz-1, 0, 0, i) = sol_part(nz-1, 0, 0, i)
308 +coef(indice)*sol_hom2(nz-1, 0, 0, i) ;
310 delete resu.set_spectral_va().c ;
311 resu.set_spectral_va().c = 0x0 ;
313 resu.set_spectral_va().ylm_i() ;
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define R_CHEB
base de Chebychev ordinaire (fin)
const Map & get_mp() const
Returns the mapping.