LORENE
get_operateur.C
1/*
2 * Copyright (c) 2000-2001 Jerome Novak
3 *
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
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 $" ;
24
25/*
26 * $Id: get_operateur.C,v 1.9 2014/10/13 08:53:28 j_novak Exp $
27 * $Log: get_operateur.C,v $
28 * Revision 1.9 2014/10/13 08:53:28 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.8 2014/10/06 15:16:08 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.7 2008/08/27 08:51:15 jl_cornou
35 * Added Jacobi(0,2) polynomials
36 *
37 * Revision 1.6 2005/01/27 10:19:43 j_novak
38 * Now using Diff operators to build the matrices.
39 *
40 * Revision 1.5 2003/06/18 08:45:27 j_novak
41 * In class Mg3d: added the member get_radial, returning only a radial grid
42 * For dAlembert solver: the way the coefficients of the operator are defined has been changed.
43 *
44 * Revision 1.4 2002/01/03 15:30:28 j_novak
45 * Some comments modified.
46 *
47 * Revision 1.3 2002/01/03 13:18:41 j_novak
48 * Optimization: the members set(i,j) and operator(i,j) of class Matrice are
49 * now defined inline. Matrice is a friend class of Tbl.
50 *
51 * Revision 1.2 2002/01/02 14:07:57 j_novak
52 * Dalembert equation is now solved in the shells. However, the number of
53 * points in theta and phi must be the same in each domain. The solver is not
54 * completely tested (beta version!).
55 *
56 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
57 * LORENE
58 *
59 * Revision 1.3 2001/10/29 10:55:28 novak
60 * Error fixed for r^2 d^2/dr^2 operator
61 *
62 * Revision 1.2 2000/12/18 13:33:46 novak
63 * *** empty log message ***
64 *
65 * Revision 1.1 2000/12/04 16:36:50 novak
66 * Initial revision
67 *
68 *
69 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/get_operateur.C,v 1.9 2014/10/13 08:53:28 j_novak Exp $
70 *
71 */
72
73// Header C :
74#include <cmath>
75
76// Headers Lorene :
77#include "param.h"
78#include "base_val.h"
79#include "diff.h"
80#include "proto.h"
81
82/*************************************************************************
83 *
84 * Routine used by sol_dalembert, to compute the matrix of the operator
85 * to be solved. The coefficients (Ci) are stored in par.get_tbl_mod(1->9)
86 * The time-step is par.get_double(0). Other inputs are:
87 * l: spherical harmonic number
88 * alpha, beta: coefficients of the affine mapping (see map.h)
89 * Outputs are: type_dal : type of the operator (see type_parite.h)
90 * operateur: matrix of the operator
91 *
92 * The operator reads:
93 *
94 * Indentity - 0.5*dt^2 [ (C1 + C3r^2) d^2/dr^2 + (C6/r + C5r) d/dr
95 * (C9/r^2 + C7) Id ]
96 *
97 *************************************************************************/
98
99
100
101 //-----------------------------------
102 // Routine pour les cas non prevus --
103 //-----------------------------------
104
105namespace Lorene {
106void _get_operateur_dal_pas_prevu(const Param& , const int&, int& , Matrice& )
107{
108 cout << "get_operateur_dal pas prevu..." << endl ;
109 abort() ;
110 exit(-1) ;
111}
112
113
114void _get_operateur_dal_r_cheb(const Param& par, const int& lz,
115int& type_dal, Matrice& operateur)
116{
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 ) ;
123
124 double dt = par.get_double(0) ;
125 dt *= 0.5*dt ;
126
127 // Copies the global coefficients to a local Tbl
128 Tbl coeff(10) ;
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) ;
141
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. ;
146
147 bool dege = (fabs(coeff(9)) < 1.e-10) ;
148 switch (dege) {
149 case true:
150 a00 = R1 - dt*(coeff(7)*R1 + coeff(8)) ;
151 a01 = R2 - dt*R2*coeff(7) ;
152 a02 = 0 ;
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) ;
156 a13 = 0 ;
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) ;
161 a24 = 0 ;
162 type_dal = ((0.1*(fabs(a20)+fabs(a21)+fabs(a22)+fabs(a23))*nr*nr*nr
163 < 1.) ? O2DEGE_SMALL : O2DEGE_LARGE ) ;
164 break ;
165 case false:
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))
179 *nr*nr*nr < 1.) ? O2NOND_SMALL : O2NOND_LARGE ) ;
180 break ;
181 }
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 ;
194
195
196
197 Diff_id id(R_CHEB, nr) ;
198 if (fabs(a00)>1.e-15) {
199 operateur = a00*id ;
200 }
201 else{
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. ;
206 }
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() ;
218
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) ;
227
228 }
229}
230
231void _get_operateur_dal_r_chebp(const Param& par, const int& lzone,
232 int& type_dal, Matrice& operateur)
233{
234 assert(lzone == 0) ; // Nucleus!
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 ) ;
241
242 double dt = par.get_double(0) ;
243 dt *= 0.5*dt ;
244
245 // Copies the global coefficients to a local Tbl and adds the -l(l+1) term
246 Tbl coeff(7) ;
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) ;
261
262 //***********************************************************************
263 // Definition of the type of operator
264 // For each type and a fixed time-step, if the number of points (nr) is too
265 // large, the round-off error makes the matrix ill-conditioned. So one has
266 // to pass the last line of the matrix to the first place (see dal_inverse).
267 // The linear combinations to put the matrix into a banded form also depend
268 // on the type of operator.
269 //***********************************************************************
270
271 if (fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(5)) < 1.e-30) {
272 // First order operator
273 if (dt < 0.1/(fabs(coeff(3)) + fabs(coeff(4))*nr))
274 type_dal = ORDRE1_SMALL ;
275 else type_dal = ORDRE1_LARGE ;
276 }
277 else {
278 // Second order degenerate (no 1/r^2 term)
279 if (fabs(coeff(5)) < 1.e-24) {
280 if (dt < 1./(fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(3))*nr
281 +fabs(coeff(4)))/nr/nr) type_dal = O2DEGE_SMALL ;
282 else type_dal = O2DEGE_LARGE ;
283 }
284 else {
285 // Second order non-degenerate (most general case)
286 if (dt < 1./((fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(3))*nr
287 + fabs(coeff(4)) + fabs(coeff(5)))*nr*nr))
288 type_dal = O2NOND_SMALL ;
289 else type_dal = O2NOND_LARGE ;
290 }
291 }
292
293 coeff.set(1) *= dt/alpha2 ;
294 coeff.set(2) *= dt ;
295 coeff.set(3) *= dt/alpha2 ;
296 coeff.set(4) *= dt ;
297 coeff.set(5) *= dt/alpha2 ;
298 coeff.set(6) *= dt ;
299
300 Diff_id id(R_CHEBP, nr) ;
301 if (fabs(1-coeff(6))>1.e-15) {
302 operateur = (1-coeff(6))*id ;
303 }
304 else{
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. ;
309 }
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() ;
315
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) ;
322 }
323
324
325}
326
327
328void _get_operateur_dal_r_chebi(const Param& par, const int& lzone,
329 int& type_dal, Matrice& operateur)
330{
331 assert(lzone == 0) ; // Nucleus!
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 ) ;
338
339 double dt = par.get_double(0) ;
340 dt *= 0.5*dt ;
341
342 // Copies the global coefficients to a local Tbl and adds the -l(l+1) term
343 Tbl coeff(7) ;
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) ;
358
359 //***********************************************************************
360 // Definition of the type of operator
361 // For each type and a fixed time-step, if the number of points (nr) is too
362 // large, the round-off error makes the matrix ill-conditioned. So one has
363 // to pass the last line of the matrix to the first place (see dal_inverse).
364 // The linear combinations to put the matrix into a banded form also depend
365 // on the type of operator.
366 //***********************************************************************
367
368 if (fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(3)) +
369 fabs(coeff(5)) < 1.e-30) {
370 // First order operator
371 if (dt < 0.1/(fabs(coeff(4))*nr))
372 type_dal = ORDRE1_SMALL ;
373 else type_dal = ORDRE1_LARGE ;
374 }
375 else {
376 if (fabs(coeff(5)+coeff(3)) < 1.e-6) {
377 // Second order degenerate (no 1/r^2 term)
378 if (dt < 1./(fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(3))*nr
379 +fabs(coeff(4)))/nr/nr) type_dal = O2DEGE_SMALL ;
380 else type_dal = O2DEGE_LARGE ;
381 }
382 else {
383 // Second order non-degenerate (most general case)
384 if (dt < 1./((fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(3))*nr
385 + fabs(coeff(4)) + fabs(coeff(5)))*nr*nr))
386 type_dal = O2NOND_SMALL ;
387 else type_dal = O2NOND_LARGE ;
388 }
389 }
390
391 coeff.set(1) *= dt/alpha2 ;
392 coeff.set(2) *= dt ;
393 coeff.set(3) *= dt/alpha2 ;
394 coeff.set(4) *= dt ;
395 coeff.set(5) *= dt/alpha2 ;
396 coeff.set(6) *= dt ;
397
398 Diff_id id(R_CHEBP, nr) ;
399 if (fabs(1-coeff(6))>1.e-15) {
400 operateur = (1-coeff(6))*id ;
401 }
402 else{
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. ;
407 }
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() ;
413
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) ;
420 }
421
422}
423
424
425
426void _get_operateur_dal_r_jaco02(const Param& par, const int& lz,
427int& type_dal, Matrice& operateur)
428{
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 ) ;
435
436 double dt = par.get_double(0) ;
437 dt *= 0.5*dt ;
438
439 // Copies the global coefficients to a local Tbl
440 Tbl coeff(10) ;
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) ;
453
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. ;
458
459 bool dege = (fabs(coeff(9)) < 1.e-10) ;
460 switch (dege) {
461 case true:
462 a00 = R1 - dt*(coeff(7)*R1 + coeff(8)) ;
463 a01 = R2 - dt*R2*coeff(7) ;
464 a02 = 0 ;
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) ;
468 a13 = 0 ;
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) ;
473 a24 = 0 ;
474 type_dal = ((0.1*(fabs(a20)+fabs(a21)+fabs(a22)+fabs(a23))*nr*nr*nr
475 < 1.) ? O2DEGE_SMALL : O2DEGE_LARGE ) ;
476 break ;
477 case false:
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))
491 *nr*nr*nr < 1.) ? O2NOND_SMALL : O2NOND_LARGE ) ;
492 break ;
493 }
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 ;
506
507
508
509 Diff_id id(R_JACO02, nr) ;
510 if (fabs(a00)>1.e-15) {
511 operateur = a00*id ;
512 }
513 else{
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. ;
518 }
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() ;
530
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) ;
539
540 }
541}
542
543
544
545
546 //--------------------------
547 //- La routine a appeler ---
548 //----------------------------
549void get_operateur_dal(const Param& par, const int& lzone,
550 const int& base_r, int& type_dal, Matrice& operateur)
551{
552
553 // Routines de derivation
554 static void (*get_operateur_dal[MAX_BASE])(const Param&,
555 const int&, int&, Matrice&) ;
556 static int nap = 0 ;
557
558 // Premier appel
559 if (nap==0) {
560 nap = 1 ;
561 for (int i=0 ; i<MAX_BASE ; i++)
562 get_operateur_dal[i] = _get_operateur_dal_pas_prevu ;
563
564 // Les routines existantes
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 ;
569 }
570
571 get_operateur_dal[base_r](par, lzone, type_dal, operateur) ;
572}
573}
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_mx2
Class for the elementary differential operator multiplication by (see the base class Diff ).
Definition diff.h:289
Lorene::Diff_mx
Class for the elementary differential operator multiplication by (see the base class Diff ).
Definition diff.h:250
Lorene::Diff_sx2
Class for the elementary differential operator division by (see the base class Diff ).
Definition diff.h:369
Lorene::Diff_sxdsdx
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:450
Lorene::Diff_x2dsdx2
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:490
Lorene::Diff_x2dsdx
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:571
Lorene::Diff_x3dsdx2
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:652
Lorene::Diff_x3dsdx
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:611
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::Matrice
Matrix handling.
Definition matrice.h:152
Lorene::Matrice::get_dim
int get_dim(int i) const
Returns the dimension of the matrix.
Definition matrice.C:260
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Tbl
Basic array class.
Definition tbl.h:161
MAX_BASE
#define MAX_BASE
Nombre max. de bases differentes.
Definition type_parite.h:144
ORDRE1_SMALL
#define ORDRE1_SMALL
Operateur du premier ordre, .
Definition type_parite.h:276
R_JACO02
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition type_parite.h:188
O2DEGE_SMALL
#define O2DEGE_SMALL
Operateur du deuxieme ordre degenere .
Definition type_parite.h:280
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition type_parite.h:170
ORDRE1_LARGE
#define ORDRE1_LARGE
Operateur du premier ordre .
Definition type_parite.h:278
O2NOND_LARGE
#define O2NOND_LARGE
Operateur du deuxieme ordre non degenere .
Definition type_parite.h:286
O2DEGE_LARGE
#define O2DEGE_LARGE
Operateur du deuxieme ordre degenere .
Definition type_parite.h:282
O2NOND_SMALL
#define O2NOND_SMALL
Operateur du deuxieme ordre non degenere .
Definition type_parite.h:284
TRA_R
#define TRA_R
Translation en R, used for a bitwise shift (in hex).
Definition type_parite.h:158
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition type_parite.h:166
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
Lorene
Lorene prototypes.
Definition app_hor.h:64