23char get_operateur_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/get_operateur.C,v 1.9 2014/10/13 08:53:28 j_novak Exp $" ;
106void _get_operateur_dal_pas_prevu(
const Param& ,
const int&,
int& , Matrice& )
108 cout <<
"get_operateur_dal pas prevu..." << endl ;
114void _get_operateur_dal_r_cheb(
const Param& par,
const int& lz,
115int& type_dal, Matrice& operateur)
117 int nr = operateur.get_dim(0) ;
118 assert (nr == operateur.get_dim(1)) ;
119 assert (par.get_n_double() > 0) ;
120 assert (par.get_n_tbl_mod() > 0) ;
121 assert ((par.get_tbl_mod()).get_dim(1) == 12 ) ;
122 assert ((par.get_tbl_mod()).get_ndim() ==2 ) ;
124 double dt = par.get_double(0) ;
129 coeff.set_etat_qcq() ;
130 coeff.set(1) = (par.get_tbl_mod())(1,lz) ;
131 coeff.set(2) = (par.get_tbl_mod())(2,lz) ;
132 coeff.set(3) = (par.get_tbl_mod())(3,lz) ;
133 coeff.set(4) = (par.get_tbl_mod())(4,lz) ;
134 coeff.set(5) = (par.get_tbl_mod())(5,lz) ;
135 coeff.set(6) = (par.get_tbl_mod())(6,lz) ;
136 coeff.set(7) = (par.get_tbl_mod())(7,lz) ;
137 coeff.set(8) = (par.get_tbl_mod())(8,lz) ;
138 coeff.set(9) = (par.get_tbl_mod())(9,lz) ;
139 double R1 = (par.get_tbl_mod())(10,lz) ;
140 double R2 = (par.get_tbl_mod())(11,lz) ;
142 double a00 = 0. ;
double a01 = 0. ;
double a02 = 0. ;
143 double a10 = 0. ;
double a11 = 0. ;
double a12 = 0. ;
144 double a13 = 0. ;
double a20 = 0. ;
double a21 = 0. ;
145 double a22 = 0. ;
double a23 = 0. ;
double a24 = 0. ;
147 bool dege = (fabs(coeff(9)) < 1.e-10) ;
150 a00 = R1 - dt*(coeff(7)*R1 + coeff(8)) ;
151 a01 = R2 - dt*R2*coeff(7) ;
153 a10 = -dt*(R1*coeff(4) + R1*R1*coeff(5) + coeff(6))/R2 ;
154 a11 = -dt*(coeff(4) + 2*R1*coeff(5)) ;
155 a12 = -dt*R2*coeff(5) ;
157 a20 = -dt*R1/(R2*R2)*(coeff(1) + R1*coeff(2) + R1*R1*coeff(3)) ;
158 a21 = -dt/R2*(coeff(1) + 2*R1*coeff(2) + 3*R1*R1*coeff(3)) ;
159 a22 = -dt*(coeff(2) + 3*R1*coeff(3)) ;
160 a23 = -dt*R2*coeff(3) ;
162 type_dal = ((0.1*(fabs(a20)+fabs(a21)+fabs(a22)+fabs(a23))*nr*nr*nr
166 a00 = R1*R1 - dt*(coeff(7)*R1*R1 + coeff(8)*R1 + coeff(9)) ;
167 a01 = 2*R1*R2 - dt*(2*R1*R2*coeff(7) + R2*coeff(8)) ;
168 a02 = R2*R2*(1 - dt*coeff(7)) ;
169 a10 = -dt*R1/R2*(R1*coeff(4) + R1*R1*coeff(5) + coeff(6)) ;
170 a11 = -dt*(2*R1*coeff(4) + 3*R1*R1*coeff(5) + coeff(6)) ;
171 a12 = -dt*(R2*coeff(4) + 3*R1*R2*coeff(5)) ;
172 a13 = -dt*R2*R2*coeff(5) ;
173 a20 = -dt*(R1*R1)/(R2*R2)*(coeff(1) + R1*coeff(2) + R1*R1*coeff(3)) ;
174 a21 = -dt*R1/R2*(2*coeff(1) + 3*R1*coeff(2) + 4*R1*R1*coeff(3)) ;
175 a22 = -dt*(coeff(1) + 3*R1*coeff(2) + 6*R1*R1*coeff(3)) ;
176 a23 = -dt*(R2*coeff(2) + 4*R1*R2*coeff(3)) ;
177 a24 = -dt*R2*R2*coeff(3) ;
178 type_dal = ((0.1*(fabs(a20)+fabs(a21)+fabs(a22)+fabs(a23)+fabs(a24))
182 if (fabs(a00)<1.e-15) a00 = 0 ;
183 if (fabs(a01)<1.e-15) a01 = 0 ;
184 if (fabs(a02)<1.e-15) a02 = 0 ;
185 if (fabs(a10)<1.e-15) a10 = 0 ;
186 if (fabs(a11)<1.e-15) a11 = 0 ;
187 if (fabs(a12)<1.e-15) a12 = 0 ;
188 if (fabs(a13)<1.e-15) a13 = 0 ;
189 if (fabs(a20)<1.e-15) a20 = 0 ;
190 if (fabs(a21)<1.e-15) a21 = 0 ;
191 if (fabs(a22)<1.e-15) a22 = 0 ;
192 if (fabs(a23)<1.e-15) a23 = 0 ;
193 if (fabs(a24)<1.e-15) a24 = 0 ;
198 if (fabs(a00)>1.e-15) {
202 operateur.set_etat_qcq() ;
203 for (
int i=0; i<nr; i++)
204 for (
int j=0; j<nr; j++)
205 operateur.set(i,j) = 0. ;
207 Diff_mx op01(
R_CHEB, nr) ;
const Matrice& m01 = op01.get_matrice() ;
208 Diff_mx2 op02(
R_CHEB, nr) ;
const Matrice& m02 = op02.get_matrice() ;
209 Diff_dsdx op10(
R_CHEB, nr) ;
const Matrice& m10 = op10.get_matrice() ;
210 Diff_xdsdx op11(
R_CHEB, nr) ;
const Matrice& m11 = op11.get_matrice() ;
211 Diff_x2dsdx op12(
R_CHEB, nr) ;
const Matrice& m12 = op12.get_matrice() ;
212 Diff_x3dsdx op13(
R_CHEB, nr) ;
const Matrice& m13 = op13.get_matrice() ;
213 Diff_dsdx2 op20(
R_CHEB, nr) ;
const Matrice& m20 = op20.get_matrice() ;
214 Diff_xdsdx2 op21(
R_CHEB, nr) ;
const Matrice& m21 = op21.get_matrice() ;
215 Diff_x2dsdx2 op22(
R_CHEB, nr) ;
const Matrice& m22 = op22.get_matrice() ;
216 Diff_x3dsdx2 op23(
R_CHEB, nr) ;
const Matrice& m23 = op23.get_matrice() ;
217 Diff_x4dsdx2 op24(
R_CHEB, nr) ;
const Matrice& m24 = op24.get_matrice() ;
219 for (
int i=0; i<nr; i++) {
220 int jmin = (i>3 ? i-3 : 0) ;
221 int jmax = (i<nr-9 ? i+10 : nr) ;
222 for (
int j=jmin ; j<jmax; j++)
223 operateur.set(i,j) += a01*m01(i,j) + a02*m02(i,j)
224 + a10*m10(i,j) + a11*m11(i,j) + a12*m12(i,j)
225 + a13*m13(i,j) + a20*m20(i,j) + a21*m21(i,j)
226 + a22*m22(i,j) + a23*m23(i,j) + a24*m24(i,j) ;
231void _get_operateur_dal_r_chebp(
const Param& par,
const int& lzone,
232 int& type_dal, Matrice& operateur)
235 int nr = operateur.get_dim(0) ;
236 assert (nr == operateur.get_dim(1)) ;
237 assert (par.get_n_double() > 0) ;
238 assert (par.get_n_tbl_mod() > 0) ;
239 assert ((par.get_tbl_mod()).get_dim(1) == 12 ) ;
240 assert ((par.get_tbl_mod()).get_ndim() ==2 ) ;
242 double dt = par.get_double(0) ;
247 coeff.set_etat_qcq() ;
248 coeff.set(1) = (par.get_tbl_mod())(1,lzone) ;
249 if (fabs(coeff(1))<1.e-15) coeff.set(1) = 0 ;
250 coeff.set(2) = (par.get_tbl_mod())(3,lzone) ;
251 if (fabs(coeff(2))<1.e-15) coeff.set(2) = 0 ;
252 coeff.set(3) = (par.get_tbl_mod())(6,lzone) ;
253 if (fabs(coeff(3))<1.e-15) coeff.set(3) = 0 ;
254 coeff.set(4) = (par.get_tbl_mod())(5,lzone) ;
255 if (fabs(coeff(4))<1.e-15) coeff.set(4) = 0 ;
256 coeff.set(5) = (par.get_tbl_mod())(9,lzone) ;
257 if (fabs(coeff(5))<1.e-15) coeff.set(5) = 0 ;
258 coeff.set(6) = (par.get_tbl_mod())(7,lzone) ;
259 if (fabs(coeff(6))<1.e-15) coeff.set(6) = 0 ;
260 double alpha2 = (par.get_tbl_mod())(11,lzone)*(par.get_tbl_mod())(11,lzone) ;
271 if (fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(5)) < 1.e-30) {
273 if (dt < 0.1/(fabs(coeff(3)) + fabs(coeff(4))*nr))
279 if (fabs(coeff(5)) < 1.e-24) {
280 if (dt < 1./(fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(3))*nr
286 if (dt < 1./((fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(3))*nr
287 + fabs(coeff(4)) + fabs(coeff(5)))*nr*nr))
293 coeff.set(1) *= dt/alpha2 ;
295 coeff.set(3) *= dt/alpha2 ;
297 coeff.set(5) *= dt/alpha2 ;
301 if (fabs(1-coeff(6))>1.e-15) {
302 operateur = (1-coeff(6))*
id ;
305 operateur.set_etat_qcq() ;
306 for (
int i=0; i<nr; i++)
307 for (
int j=0; j<nr; j++)
308 operateur.set(i,j) = 0. ;
310 Diff_sx2 op02(
R_CHEBP, nr) ;
const Matrice& m02 = op02.get_matrice() ;
311 Diff_xdsdx op11(
R_CHEBP, nr) ;
const Matrice& m11 = op11.get_matrice() ;
312 Diff_sxdsdx op12(
R_CHEBP, nr) ;
const Matrice& m12 = op12.get_matrice() ;
313 Diff_dsdx2 op20(
R_CHEBP, nr) ;
const Matrice& m20 = op20.get_matrice() ;
314 Diff_x2dsdx2 op22(
R_CHEBP, nr) ;
const Matrice& m22 = op22.get_matrice() ;
316 for (
int i=0; i<nr; i++) {
317 int jmin = (i>3 ? i-3 : 0) ;
318 int jmax = (i<nr-9 ? i+10 : nr) ;
319 for (
int j=jmin ; j<jmax; j++)
320 operateur.set(i,j) -= coeff(1)*m20(i,j) + coeff(2)*m22(i,j)
321 + coeff(3)*m12(i,j) + coeff(4)*m11(i,j) + coeff(5)*m02(i,j) ;
328void _get_operateur_dal_r_chebi(
const Param& par,
const int& lzone,
329 int& type_dal, Matrice& operateur)
332 int nr = operateur.get_dim(0) ;
333 assert (nr == operateur.get_dim(1)) ;
334 assert (par.get_n_double() > 0) ;
335 assert (par.get_n_tbl_mod() > 0) ;
336 assert ((par.get_tbl_mod()).get_dim(1) == 12 ) ;
337 assert ((par.get_tbl_mod()).get_ndim() == 2 ) ;
339 double dt = par.get_double(0) ;
344 coeff.set_etat_qcq() ;
345 coeff.set(1) = (par.get_tbl_mod())(1,lzone) ;
346 if (fabs(coeff(1))<1.e-15) coeff.set(1) = 0 ;
347 coeff.set(2) = (par.get_tbl_mod())(3,lzone) ;
348 if (fabs(coeff(2))<1.e-15) coeff.set(2) = 0 ;
349 coeff.set(3) = (par.get_tbl_mod())(6,lzone) ;
350 if (fabs(coeff(3))<1.e-15) coeff.set(3) = 0 ;
351 coeff.set(4) = (par.get_tbl_mod())(5,lzone) ;
352 if (fabs(coeff(4))<1.e-15) coeff.set(4) = 0 ;
353 coeff.set(5) = (par.get_tbl_mod())(9,lzone) ;
354 if (fabs(coeff(5))<1.e-15) coeff.set(5) = 0 ;
355 coeff.set(6) = (par.get_tbl_mod())(7,lzone) ;
356 if (fabs(coeff(6))<1.e-15) coeff.set(6) = 0 ;
357 double alpha2 = (par.get_tbl_mod())(11,lzone)*(par.get_tbl_mod())(11,lzone) ;
368 if (fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(3)) +
369 fabs(coeff(5)) < 1.e-30) {
371 if (dt < 0.1/(fabs(coeff(4))*nr))
376 if (fabs(coeff(5)+coeff(3)) < 1.e-6) {
378 if (dt < 1./(fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(3))*nr
384 if (dt < 1./((fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(3))*nr
385 + fabs(coeff(4)) + fabs(coeff(5)))*nr*nr))
391 coeff.set(1) *= dt/alpha2 ;
393 coeff.set(3) *= dt/alpha2 ;
395 coeff.set(5) *= dt/alpha2 ;
399 if (fabs(1-coeff(6))>1.e-15) {
400 operateur = (1-coeff(6))*
id ;
403 operateur.set_etat_qcq() ;
404 for (
int i=0; i<nr; i++)
405 for (
int j=0; j<nr; j++)
406 operateur.set(i,j) = 0. ;
408 Diff_sx2 op02(
R_CHEBI, nr) ;
const Matrice& m02 = op02.get_matrice() ;
409 Diff_xdsdx op11(
R_CHEBI, nr) ;
const Matrice& m11 = op11.get_matrice() ;
410 Diff_sxdsdx op12(
R_CHEBI, nr) ;
const Matrice& m12 = op12.get_matrice() ;
411 Diff_dsdx2 op20(
R_CHEBI, nr) ;
const Matrice& m20 = op20.get_matrice() ;
412 Diff_x2dsdx2 op22(
R_CHEBI, nr) ;
const Matrice& m22 = op22.get_matrice() ;
414 for (
int i=0; i<nr; i++) {
415 int jmin = (i>3 ? i-3 : 0) ;
416 int jmax = (i<nr-9 ? i+10 : nr) ;
417 for (
int j=jmin ; j<jmax; j++)
418 operateur.set(i,j) -= coeff(1)*m20(i,j) + coeff(2)*m22(i,j)
419 + coeff(3)*m12(i,j) + coeff(4)*m11(i,j) + coeff(5)*m02(i,j) ;
426void _get_operateur_dal_r_jaco02(
const Param& par,
const int& lz,
427int& type_dal, Matrice& operateur)
429 int nr = operateur.get_dim(0) ;
430 assert (nr == operateur.get_dim(1)) ;
431 assert (par.get_n_double() > 0) ;
432 assert (par.get_n_tbl_mod() > 0) ;
433 assert ((par.get_tbl_mod()).get_dim(1) == 12 ) ;
434 assert ((par.get_tbl_mod()).get_ndim() ==2 ) ;
436 double dt = par.get_double(0) ;
441 coeff.set_etat_qcq() ;
442 coeff.set(1) = (par.get_tbl_mod())(1,lz) ;
443 coeff.set(2) = (par.get_tbl_mod())(2,lz) ;
444 coeff.set(3) = (par.get_tbl_mod())(3,lz) ;
445 coeff.set(4) = (par.get_tbl_mod())(4,lz) ;
446 coeff.set(5) = (par.get_tbl_mod())(5,lz) ;
447 coeff.set(6) = (par.get_tbl_mod())(6,lz) ;
448 coeff.set(7) = (par.get_tbl_mod())(7,lz) ;
449 coeff.set(8) = (par.get_tbl_mod())(8,lz) ;
450 coeff.set(9) = (par.get_tbl_mod())(9,lz) ;
451 double R1 = (par.get_tbl_mod())(10,lz) ;
452 double R2 = (par.get_tbl_mod())(11,lz) ;
454 double a00 = 0. ;
double a01 = 0. ;
double a02 = 0. ;
455 double a10 = 0. ;
double a11 = 0. ;
double a12 = 0. ;
456 double a13 = 0. ;
double a20 = 0. ;
double a21 = 0. ;
457 double a22 = 0. ;
double a23 = 0. ;
double a24 = 0. ;
459 bool dege = (fabs(coeff(9)) < 1.e-10) ;
462 a00 = R1 - dt*(coeff(7)*R1 + coeff(8)) ;
463 a01 = R2 - dt*R2*coeff(7) ;
465 a10 = -dt*(R1*coeff(4) + R1*R1*coeff(5) + coeff(6))/R2 ;
466 a11 = -dt*(coeff(4) + 2*R1*coeff(5)) ;
467 a12 = -dt*R2*coeff(5) ;
469 a20 = -dt*R1/(R2*R2)*(coeff(1) + R1*coeff(2) + R1*R1*coeff(3)) ;
470 a21 = -dt/R2*(coeff(1) + 2*R1*coeff(2) + 3*R1*R1*coeff(3)) ;
471 a22 = -dt*(coeff(2) + 3*R1*coeff(3)) ;
472 a23 = -dt*R2*coeff(3) ;
474 type_dal = ((0.1*(fabs(a20)+fabs(a21)+fabs(a22)+fabs(a23))*nr*nr*nr
478 a00 = R1*R1 - dt*(coeff(7)*R1*R1 + coeff(8)*R1 + coeff(9)) ;
479 a01 = 2*R1*R2 - dt*(2*R1*R2*coeff(7) + R2*coeff(8)) ;
480 a02 = R2*R2*(1 - dt*coeff(7)) ;
481 a10 = -dt*R1/R2*(R1*coeff(4) + R1*R1*coeff(5) + coeff(6)) ;
482 a11 = -dt*(2*R1*coeff(4) + 3*R1*R1*coeff(5) + coeff(6)) ;
483 a12 = -dt*(R2*coeff(4) + 3*R1*R2*coeff(5)) ;
484 a13 = -dt*R2*R2*coeff(5) ;
485 a20 = -dt*(R1*R1)/(R2*R2)*(coeff(1) + R1*coeff(2) + R1*R1*coeff(3)) ;
486 a21 = -dt*R1/R2*(2*coeff(1) + 3*R1*coeff(2) + 4*R1*R1*coeff(3)) ;
487 a22 = -dt*(coeff(1) + 3*R1*coeff(2) + 6*R1*R1*coeff(3)) ;
488 a23 = -dt*(R2*coeff(2) + 4*R1*R2*coeff(3)) ;
489 a24 = -dt*R2*R2*coeff(3) ;
490 type_dal = ((0.1*(fabs(a20)+fabs(a21)+fabs(a22)+fabs(a23)+fabs(a24))
494 if (fabs(a00)<1.e-15) a00 = 0 ;
495 if (fabs(a01)<1.e-15) a01 = 0 ;
496 if (fabs(a02)<1.e-15) a02 = 0 ;
497 if (fabs(a10)<1.e-15) a10 = 0 ;
498 if (fabs(a11)<1.e-15) a11 = 0 ;
499 if (fabs(a12)<1.e-15) a12 = 0 ;
500 if (fabs(a13)<1.e-15) a13 = 0 ;
501 if (fabs(a20)<1.e-15) a20 = 0 ;
502 if (fabs(a21)<1.e-15) a21 = 0 ;
503 if (fabs(a22)<1.e-15) a22 = 0 ;
504 if (fabs(a23)<1.e-15) a23 = 0 ;
505 if (fabs(a24)<1.e-15) a24 = 0 ;
510 if (fabs(a00)>1.e-15) {
514 operateur.set_etat_qcq() ;
515 for (
int i=0; i<nr; i++)
516 for (
int j=0; j<nr; j++)
517 operateur.set(i,j) = 0. ;
519 Diff_mx op01(
R_JACO02, nr) ;
const Matrice& m01 = op01.get_matrice() ;
520 Diff_mx2 op02(
R_JACO02, nr) ;
const Matrice& m02 = op02.get_matrice() ;
521 Diff_dsdx op10(
R_JACO02, nr) ;
const Matrice& m10 = op10.get_matrice() ;
522 Diff_xdsdx op11(
R_JACO02, nr) ;
const Matrice& m11 = op11.get_matrice() ;
523 Diff_x2dsdx op12(
R_JACO02, nr) ;
const Matrice& m12 = op12.get_matrice() ;
524 Diff_x3dsdx op13(
R_JACO02, nr) ;
const Matrice& m13 = op13.get_matrice() ;
525 Diff_dsdx2 op20(
R_JACO02, nr) ;
const Matrice& m20 = op20.get_matrice() ;
526 Diff_xdsdx2 op21(
R_JACO02, nr) ;
const Matrice& m21 = op21.get_matrice() ;
527 Diff_x2dsdx2 op22(
R_JACO02, nr) ;
const Matrice& m22 = op22.get_matrice() ;
528 Diff_x3dsdx2 op23(
R_JACO02, nr) ;
const Matrice& m23 = op23.get_matrice() ;
529 Diff_x4dsdx2 op24(
R_JACO02, nr) ;
const Matrice& m24 = op24.get_matrice() ;
531 for (
int i=0; i<nr; i++) {
532 int jmin = (i>3 ? i-3 : 0) ;
533 int jmax = (i<nr-9 ? i+10 : nr) ;
534 for (
int j=jmin ; j<jmax; j++)
535 operateur.set(i,j) += a01*m01(i,j) + a02*m02(i,j)
536 + a10*m10(i,j) + a11*m11(i,j) + a12*m12(i,j)
537 + a13*m13(i,j) + a20*m20(i,j) + a21*m21(i,j)
538 + a22*m22(i,j) + a23*m23(i,j) + a24*m24(i,j) ;
549void get_operateur_dal(
const Param& par,
const int& lzone,
550 const int& base_r,
int& type_dal, Matrice& operateur)
554 static void (*get_operateur_dal[
MAX_BASE])(
const Param&,
555 const int&,
int&, Matrice&) ;
562 get_operateur_dal[i] = _get_operateur_dal_pas_prevu ;
565 get_operateur_dal[
R_CHEB >>
TRA_R] = _get_operateur_dal_r_cheb ;
566 get_operateur_dal[
R_CHEBP >>
TRA_R] = _get_operateur_dal_r_chebp ;
567 get_operateur_dal[
R_CHEBI >>
TRA_R] = _get_operateur_dal_r_chebi ;
568 get_operateur_dal[
R_JACO02 >>
TRA_R] = _get_operateur_dal_r_jaco02 ;
571 get_operateur_dal[base_r](par, lzone, type_dal, operateur) ;
#define MAX_BASE
Nombre max. de bases differentes.
#define ORDRE1_SMALL
Operateur du premier ordre, .
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
#define O2DEGE_SMALL
Operateur du deuxieme ordre degenere .
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define ORDRE1_LARGE
Operateur du premier ordre .
#define O2NOND_LARGE
Operateur du deuxieme ordre non degenere .
#define O2DEGE_LARGE
Operateur du deuxieme ordre degenere .
#define O2NOND_SMALL
Operateur du deuxieme ordre non degenere .
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
#define R_CHEB
base de Chebychev ordinaire (fin)
#define R_CHEBP
base de Cheb. paire (rare) seulement