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
23
24
25/*
26 * $Id: get_operateur.C,v 1.10 2016/12/05 16:18:09 j_novak Exp $
27 * $Log: get_operateur.C,v $
28 * Revision 1.10 2016/12/05 16:18:09 j_novak
29 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
30 *
31 * Revision 1.9 2014/10/13 08:53:28 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.8 2014/10/06 15:16:08 j_novak
35 * Modified #include directives to use c++ syntax.
36 *
37 * Revision 1.7 2008/08/27 08:51:15 jl_cornou
38 * Added Jacobi(0,2) polynomials
39 *
40 * Revision 1.6 2005/01/27 10:19:43 j_novak
41 * Now using Diff operators to build the matrices.
42 *
43 * Revision 1.5 2003/06/18 08:45:27 j_novak
44 * In class Mg3d: added the member get_radial, returning only a radial grid
45 * For dAlembert solver: the way the coefficients of the operator are defined has been changed.
46 *
47 * Revision 1.4 2002/01/03 15:30:28 j_novak
48 * Some comments modified.
49 *
50 * Revision 1.3 2002/01/03 13:18:41 j_novak
51 * Optimization: the members set(i,j) and operator(i,j) of class Matrice are
52 * now defined inline. Matrice is a friend class of Tbl.
53 *
54 * Revision 1.2 2002/01/02 14:07:57 j_novak
55 * Dalembert equation is now solved in the shells. However, the number of
56 * points in theta and phi must be the same in each domain. The solver is not
57 * completely tested (beta version!).
58 *
59 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
60 * LORENE
61 *
62 * Revision 1.3 2001/10/29 10:55:28 novak
63 * Error fixed for r^2 d^2/dr^2 operator
64 *
65 * Revision 1.2 2000/12/18 13:33:46 novak
66 * *** empty log message ***
67 *
68 * Revision 1.1 2000/12/04 16:36:50 novak
69 * Initial revision
70 *
71 *
72 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/get_operateur.C,v 1.10 2016/12/05 16:18:09 j_novak Exp $
73 *
74 */
75
76// Header C :
77#include <cmath>
78
79// Headers Lorene :
80#include "param.h"
81#include "base_val.h"
82#include "diff.h"
83#include "proto.h"
84
85/*************************************************************************
86 *
87 * Routine used by sol_dalembert, to compute the matrix of the operator
88 * to be solved. The coefficients (Ci) are stored in par.get_tbl_mod(1->9)
89 * The time-step is par.get_double(0). Other inputs are:
90 * l: spherical harmonic number
91 * alpha, beta: coefficients of the affine mapping (see map.h)
92 * Outputs are: type_dal : type of the operator (see type_parite.h)
93 * operateur: matrix of the operator
94 *
95 * The operator reads:
96 *
97 * Indentity - 0.5*dt^2 [ (C1 + C3r^2) d^2/dr^2 + (C6/r + C5r) d/dr
98 * (C9/r^2 + C7) Id ]
99 *
100 *************************************************************************/
101
102
103
104 //-----------------------------------
105 // Routine pour les cas non prevus --
106 //-----------------------------------
107
108namespace Lorene {
109void _get_operateur_dal_pas_prevu(const Param& , const int&, int& , Matrice& )
110{
111 cout << "get_operateur_dal pas prevu..." << endl ;
112 abort() ;
113 exit(-1) ;
114}
115
116
117void _get_operateur_dal_r_cheb(const Param& par, const int& lz,
118int& type_dal, Matrice& operateur)
119{
120 int nr = operateur.get_dim(0) ;
121 assert (nr == operateur.get_dim(1)) ;
122 assert (par.get_n_double() > 0) ;
123 assert (par.get_n_tbl_mod() > 0) ;
124 assert ((par.get_tbl_mod()).get_dim(1) == 12 ) ;
125 assert ((par.get_tbl_mod()).get_ndim() ==2 ) ;
126
127 double dt = par.get_double(0) ;
128 dt *= 0.5*dt ;
129
130 // Copies the global coefficients to a local Tbl
131 Tbl coeff(10) ;
132 coeff.set_etat_qcq() ;
133 coeff.set(1) = (par.get_tbl_mod())(1,lz) ;
134 coeff.set(2) = (par.get_tbl_mod())(2,lz) ;
135 coeff.set(3) = (par.get_tbl_mod())(3,lz) ;
136 coeff.set(4) = (par.get_tbl_mod())(4,lz) ;
137 coeff.set(5) = (par.get_tbl_mod())(5,lz) ;
138 coeff.set(6) = (par.get_tbl_mod())(6,lz) ;
139 coeff.set(7) = (par.get_tbl_mod())(7,lz) ;
140 coeff.set(8) = (par.get_tbl_mod())(8,lz) ;
141 coeff.set(9) = (par.get_tbl_mod())(9,lz) ;
142 double R1 = (par.get_tbl_mod())(10,lz) ;
143 double R2 = (par.get_tbl_mod())(11,lz) ;
144
145 double a00 = 0. ; double a01 = 0. ; double a02 = 0. ;
146 double a10 = 0. ; double a11 = 0. ; double a12 = 0. ;
147 double a13 = 0. ; double a20 = 0. ; double a21 = 0. ;
148 double a22 = 0. ; double a23 = 0. ; double a24 = 0. ;
149
150 bool dege = (fabs(coeff(9)) < 1.e-10) ;
151 switch (dege) {
152 case true:
153 a00 = R1 - dt*(coeff(7)*R1 + coeff(8)) ;
154 a01 = R2 - dt*R2*coeff(7) ;
155 a02 = 0 ;
156 a10 = -dt*(R1*coeff(4) + R1*R1*coeff(5) + coeff(6))/R2 ;
157 a11 = -dt*(coeff(4) + 2*R1*coeff(5)) ;
158 a12 = -dt*R2*coeff(5) ;
159 a13 = 0 ;
160 a20 = -dt*R1/(R2*R2)*(coeff(1) + R1*coeff(2) + R1*R1*coeff(3)) ;
161 a21 = -dt/R2*(coeff(1) + 2*R1*coeff(2) + 3*R1*R1*coeff(3)) ;
162 a22 = -dt*(coeff(2) + 3*R1*coeff(3)) ;
163 a23 = -dt*R2*coeff(3) ;
164 a24 = 0 ;
165 type_dal = ((0.1*(fabs(a20)+fabs(a21)+fabs(a22)+fabs(a23))*nr*nr*nr
166 < 1.) ? O2DEGE_SMALL : O2DEGE_LARGE ) ;
167 break ;
168 case false:
169 a00 = R1*R1 - dt*(coeff(7)*R1*R1 + coeff(8)*R1 + coeff(9)) ;
170 a01 = 2*R1*R2 - dt*(2*R1*R2*coeff(7) + R2*coeff(8)) ;
171 a02 = R2*R2*(1 - dt*coeff(7)) ;
172 a10 = -dt*R1/R2*(R1*coeff(4) + R1*R1*coeff(5) + coeff(6)) ;
173 a11 = -dt*(2*R1*coeff(4) + 3*R1*R1*coeff(5) + coeff(6)) ;
174 a12 = -dt*(R2*coeff(4) + 3*R1*R2*coeff(5)) ;
175 a13 = -dt*R2*R2*coeff(5) ;
176 a20 = -dt*(R1*R1)/(R2*R2)*(coeff(1) + R1*coeff(2) + R1*R1*coeff(3)) ;
177 a21 = -dt*R1/R2*(2*coeff(1) + 3*R1*coeff(2) + 4*R1*R1*coeff(3)) ;
178 a22 = -dt*(coeff(1) + 3*R1*coeff(2) + 6*R1*R1*coeff(3)) ;
179 a23 = -dt*(R2*coeff(2) + 4*R1*R2*coeff(3)) ;
180 a24 = -dt*R2*R2*coeff(3) ;
181 type_dal = ((0.1*(fabs(a20)+fabs(a21)+fabs(a22)+fabs(a23)+fabs(a24))
182 *nr*nr*nr < 1.) ? O2NOND_SMALL : O2NOND_LARGE ) ;
183 break ;
184 }
185 if (fabs(a00)<1.e-15) a00 = 0 ;
186 if (fabs(a01)<1.e-15) a01 = 0 ;
187 if (fabs(a02)<1.e-15) a02 = 0 ;
188 if (fabs(a10)<1.e-15) a10 = 0 ;
189 if (fabs(a11)<1.e-15) a11 = 0 ;
190 if (fabs(a12)<1.e-15) a12 = 0 ;
191 if (fabs(a13)<1.e-15) a13 = 0 ;
192 if (fabs(a20)<1.e-15) a20 = 0 ;
193 if (fabs(a21)<1.e-15) a21 = 0 ;
194 if (fabs(a22)<1.e-15) a22 = 0 ;
195 if (fabs(a23)<1.e-15) a23 = 0 ;
196 if (fabs(a24)<1.e-15) a24 = 0 ;
197
198
199
200 Diff_id id(R_CHEB, nr) ;
201 if (fabs(a00)>1.e-15) {
202 operateur = a00*id ;
203 }
204 else{
205 operateur.set_etat_qcq() ;
206 for (int i=0; i<nr; i++)
207 for (int j=0; j<nr; j++)
208 operateur.set(i,j) = 0. ;
209 }
210 Diff_mx op01(R_CHEB, nr) ; const Matrice& m01 = op01.get_matrice() ;
211 Diff_mx2 op02(R_CHEB, nr) ; const Matrice& m02 = op02.get_matrice() ;
212 Diff_dsdx op10(R_CHEB, nr) ; const Matrice& m10 = op10.get_matrice() ;
213 Diff_xdsdx op11(R_CHEB, nr) ; const Matrice& m11 = op11.get_matrice() ;
214 Diff_x2dsdx op12(R_CHEB, nr) ; const Matrice& m12 = op12.get_matrice() ;
215 Diff_x3dsdx op13(R_CHEB, nr) ; const Matrice& m13 = op13.get_matrice() ;
216 Diff_dsdx2 op20(R_CHEB, nr) ; const Matrice& m20 = op20.get_matrice() ;
217 Diff_xdsdx2 op21(R_CHEB, nr) ; const Matrice& m21 = op21.get_matrice() ;
218 Diff_x2dsdx2 op22(R_CHEB, nr) ; const Matrice& m22 = op22.get_matrice() ;
219 Diff_x3dsdx2 op23(R_CHEB, nr) ; const Matrice& m23 = op23.get_matrice() ;
220 Diff_x4dsdx2 op24(R_CHEB, nr) ; const Matrice& m24 = op24.get_matrice() ;
221
222 for (int i=0; i<nr; i++) {
223 int jmin = (i>3 ? i-3 : 0) ;
224 int jmax = (i<nr-9 ? i+10 : nr) ;
225 for (int j=jmin ; j<jmax; j++)
226 operateur.set(i,j) += a01*m01(i,j) + a02*m02(i,j)
227 + a10*m10(i,j) + a11*m11(i,j) + a12*m12(i,j)
228 + a13*m13(i,j) + a20*m20(i,j) + a21*m21(i,j)
229 + a22*m22(i,j) + a23*m23(i,j) + a24*m24(i,j) ;
230
231 }
232}
233
234void _get_operateur_dal_r_chebp(const Param& par, const int& lzone,
235 int& type_dal, Matrice& operateur)
236{
237 assert(lzone == 0) ; // Nucleus!
238 int nr = operateur.get_dim(0) ;
239 assert (nr == operateur.get_dim(1)) ;
240 assert (par.get_n_double() > 0) ;
241 assert (par.get_n_tbl_mod() > 0) ;
242 assert ((par.get_tbl_mod()).get_dim(1) == 12 ) ;
243 assert ((par.get_tbl_mod()).get_ndim() ==2 ) ;
244
245 double dt = par.get_double(0) ;
246 dt *= 0.5*dt ;
247
248 // Copies the global coefficients to a local Tbl and adds the -l(l+1) term
249 Tbl coeff(7) ;
250 coeff.set_etat_qcq() ;
251 coeff.set(1) = (par.get_tbl_mod())(1,lzone) ;
252 if (fabs(coeff(1))<1.e-15) coeff.set(1) = 0 ;
253 coeff.set(2) = (par.get_tbl_mod())(3,lzone) ;
254 if (fabs(coeff(2))<1.e-15) coeff.set(2) = 0 ;
255 coeff.set(3) = (par.get_tbl_mod())(6,lzone) ;
256 if (fabs(coeff(3))<1.e-15) coeff.set(3) = 0 ;
257 coeff.set(4) = (par.get_tbl_mod())(5,lzone) ;
258 if (fabs(coeff(4))<1.e-15) coeff.set(4) = 0 ;
259 coeff.set(5) = (par.get_tbl_mod())(9,lzone) ;
260 if (fabs(coeff(5))<1.e-15) coeff.set(5) = 0 ;
261 coeff.set(6) = (par.get_tbl_mod())(7,lzone) ;
262 if (fabs(coeff(6))<1.e-15) coeff.set(6) = 0 ;
263 double alpha2 = (par.get_tbl_mod())(11,lzone)*(par.get_tbl_mod())(11,lzone) ;
264
265 //***********************************************************************
266 // Definition of the type of operator
267 // For each type and a fixed time-step, if the number of points (nr) is too
268 // large, the round-off error makes the matrix ill-conditioned. So one has
269 // to pass the last line of the matrix to the first place (see dal_inverse).
270 // The linear combinations to put the matrix into a banded form also depend
271 // on the type of operator.
272 //***********************************************************************
273
274 if (fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(5)) < 1.e-30) {
275 // First order operator
276 if (dt < 0.1/(fabs(coeff(3)) + fabs(coeff(4))*nr))
277 type_dal = ORDRE1_SMALL ;
278 else type_dal = ORDRE1_LARGE ;
279 }
280 else {
281 // Second order degenerate (no 1/r^2 term)
282 if (fabs(coeff(5)) < 1.e-24) {
283 if (dt < 1./(fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(3))*nr
284 +fabs(coeff(4)))/nr/nr) type_dal = O2DEGE_SMALL ;
285 else type_dal = O2DEGE_LARGE ;
286 }
287 else {
288 // Second order non-degenerate (most general case)
289 if (dt < 1./((fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(3))*nr
290 + fabs(coeff(4)) + fabs(coeff(5)))*nr*nr))
291 type_dal = O2NOND_SMALL ;
292 else type_dal = O2NOND_LARGE ;
293 }
294 }
295
296 coeff.set(1) *= dt/alpha2 ;
297 coeff.set(2) *= dt ;
298 coeff.set(3) *= dt/alpha2 ;
299 coeff.set(4) *= dt ;
300 coeff.set(5) *= dt/alpha2 ;
301 coeff.set(6) *= dt ;
302
303 Diff_id id(R_CHEBP, nr) ;
304 if (fabs(1-coeff(6))>1.e-15) {
305 operateur = (1-coeff(6))*id ;
306 }
307 else{
308 operateur.set_etat_qcq() ;
309 for (int i=0; i<nr; i++)
310 for (int j=0; j<nr; j++)
311 operateur.set(i,j) = 0. ;
312 }
313 Diff_sx2 op02(R_CHEBP, nr) ; const Matrice& m02 = op02.get_matrice() ;
314 Diff_xdsdx op11(R_CHEBP, nr) ; const Matrice& m11 = op11.get_matrice() ;
315 Diff_sxdsdx op12(R_CHEBP, nr) ; const Matrice& m12 = op12.get_matrice() ;
316 Diff_dsdx2 op20(R_CHEBP, nr) ; const Matrice& m20 = op20.get_matrice() ;
317 Diff_x2dsdx2 op22(R_CHEBP, nr) ; const Matrice& m22 = op22.get_matrice() ;
318
319 for (int i=0; i<nr; i++) {
320 int jmin = (i>3 ? i-3 : 0) ;
321 int jmax = (i<nr-9 ? i+10 : nr) ;
322 for (int j=jmin ; j<jmax; j++)
323 operateur.set(i,j) -= coeff(1)*m20(i,j) + coeff(2)*m22(i,j)
324 + coeff(3)*m12(i,j) + coeff(4)*m11(i,j) + coeff(5)*m02(i,j) ;
325 }
326
327
328}
329
330
331void _get_operateur_dal_r_chebi(const Param& par, const int& lzone,
332 int& type_dal, Matrice& operateur)
333{
334 assert(lzone == 0) ; // Nucleus!
335 int nr = operateur.get_dim(0) ;
336 assert (nr == operateur.get_dim(1)) ;
337 assert (par.get_n_double() > 0) ;
338 assert (par.get_n_tbl_mod() > 0) ;
339 assert ((par.get_tbl_mod()).get_dim(1) == 12 ) ;
340 assert ((par.get_tbl_mod()).get_ndim() == 2 ) ;
341
342 double dt = par.get_double(0) ;
343 dt *= 0.5*dt ;
344
345 // Copies the global coefficients to a local Tbl and adds the -l(l+1) term
346 Tbl coeff(7) ;
347 coeff.set_etat_qcq() ;
348 coeff.set(1) = (par.get_tbl_mod())(1,lzone) ;
349 if (fabs(coeff(1))<1.e-15) coeff.set(1) = 0 ;
350 coeff.set(2) = (par.get_tbl_mod())(3,lzone) ;
351 if (fabs(coeff(2))<1.e-15) coeff.set(2) = 0 ;
352 coeff.set(3) = (par.get_tbl_mod())(6,lzone) ;
353 if (fabs(coeff(3))<1.e-15) coeff.set(3) = 0 ;
354 coeff.set(4) = (par.get_tbl_mod())(5,lzone) ;
355 if (fabs(coeff(4))<1.e-15) coeff.set(4) = 0 ;
356 coeff.set(5) = (par.get_tbl_mod())(9,lzone) ;
357 if (fabs(coeff(5))<1.e-15) coeff.set(5) = 0 ;
358 coeff.set(6) = (par.get_tbl_mod())(7,lzone) ;
359 if (fabs(coeff(6))<1.e-15) coeff.set(6) = 0 ;
360 double alpha2 = (par.get_tbl_mod())(11,lzone)*(par.get_tbl_mod())(11,lzone) ;
361
362 //***********************************************************************
363 // Definition of the type of operator
364 // For each type and a fixed time-step, if the number of points (nr) is too
365 // large, the round-off error makes the matrix ill-conditioned. So one has
366 // to pass the last line of the matrix to the first place (see dal_inverse).
367 // The linear combinations to put the matrix into a banded form also depend
368 // on the type of operator.
369 //***********************************************************************
370
371 if (fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(3)) +
372 fabs(coeff(5)) < 1.e-30) {
373 // First order operator
374 if (dt < 0.1/(fabs(coeff(4))*nr))
375 type_dal = ORDRE1_SMALL ;
376 else type_dal = ORDRE1_LARGE ;
377 }
378 else {
379 if (fabs(coeff(5)+coeff(3)) < 1.e-6) {
380 // Second order degenerate (no 1/r^2 term)
381 if (dt < 1./(fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(3))*nr
382 +fabs(coeff(4)))/nr/nr) type_dal = O2DEGE_SMALL ;
383 else type_dal = O2DEGE_LARGE ;
384 }
385 else {
386 // Second order non-degenerate (most general case)
387 if (dt < 1./((fabs(coeff(1)) + fabs(coeff(2)) + fabs(coeff(3))*nr
388 + fabs(coeff(4)) + fabs(coeff(5)))*nr*nr))
389 type_dal = O2NOND_SMALL ;
390 else type_dal = O2NOND_LARGE ;
391 }
392 }
393
394 coeff.set(1) *= dt/alpha2 ;
395 coeff.set(2) *= dt ;
396 coeff.set(3) *= dt/alpha2 ;
397 coeff.set(4) *= dt ;
398 coeff.set(5) *= dt/alpha2 ;
399 coeff.set(6) *= dt ;
400
401 Diff_id id(R_CHEBP, nr) ;
402 if (fabs(1-coeff(6))>1.e-15) {
403 operateur = (1-coeff(6))*id ;
404 }
405 else{
406 operateur.set_etat_qcq() ;
407 for (int i=0; i<nr; i++)
408 for (int j=0; j<nr; j++)
409 operateur.set(i,j) = 0. ;
410 }
411 Diff_sx2 op02(R_CHEBI, nr) ; const Matrice& m02 = op02.get_matrice() ;
412 Diff_xdsdx op11(R_CHEBI, nr) ; const Matrice& m11 = op11.get_matrice() ;
413 Diff_sxdsdx op12(R_CHEBI, nr) ; const Matrice& m12 = op12.get_matrice() ;
414 Diff_dsdx2 op20(R_CHEBI, nr) ; const Matrice& m20 = op20.get_matrice() ;
415 Diff_x2dsdx2 op22(R_CHEBI, nr) ; const Matrice& m22 = op22.get_matrice() ;
416
417 for (int i=0; i<nr; i++) {
418 int jmin = (i>3 ? i-3 : 0) ;
419 int jmax = (i<nr-9 ? i+10 : nr) ;
420 for (int j=jmin ; j<jmax; j++)
421 operateur.set(i,j) -= coeff(1)*m20(i,j) + coeff(2)*m22(i,j)
422 + coeff(3)*m12(i,j) + coeff(4)*m11(i,j) + coeff(5)*m02(i,j) ;
423 }
424
425}
426
427
428
429void _get_operateur_dal_r_jaco02(const Param& par, const int& lz,
430int& type_dal, Matrice& operateur)
431{
432 int nr = operateur.get_dim(0) ;
433 assert (nr == operateur.get_dim(1)) ;
434 assert (par.get_n_double() > 0) ;
435 assert (par.get_n_tbl_mod() > 0) ;
436 assert ((par.get_tbl_mod()).get_dim(1) == 12 ) ;
437 assert ((par.get_tbl_mod()).get_ndim() ==2 ) ;
438
439 double dt = par.get_double(0) ;
440 dt *= 0.5*dt ;
441
442 // Copies the global coefficients to a local Tbl
443 Tbl coeff(10) ;
444 coeff.set_etat_qcq() ;
445 coeff.set(1) = (par.get_tbl_mod())(1,lz) ;
446 coeff.set(2) = (par.get_tbl_mod())(2,lz) ;
447 coeff.set(3) = (par.get_tbl_mod())(3,lz) ;
448 coeff.set(4) = (par.get_tbl_mod())(4,lz) ;
449 coeff.set(5) = (par.get_tbl_mod())(5,lz) ;
450 coeff.set(6) = (par.get_tbl_mod())(6,lz) ;
451 coeff.set(7) = (par.get_tbl_mod())(7,lz) ;
452 coeff.set(8) = (par.get_tbl_mod())(8,lz) ;
453 coeff.set(9) = (par.get_tbl_mod())(9,lz) ;
454 double R1 = (par.get_tbl_mod())(10,lz) ;
455 double R2 = (par.get_tbl_mod())(11,lz) ;
456
457 double a00 = 0. ; double a01 = 0. ; double a02 = 0. ;
458 double a10 = 0. ; double a11 = 0. ; double a12 = 0. ;
459 double a13 = 0. ; double a20 = 0. ; double a21 = 0. ;
460 double a22 = 0. ; double a23 = 0. ; double a24 = 0. ;
461
462 bool dege = (fabs(coeff(9)) < 1.e-10) ;
463 switch (dege) {
464 case true:
465 a00 = R1 - dt*(coeff(7)*R1 + coeff(8)) ;
466 a01 = R2 - dt*R2*coeff(7) ;
467 a02 = 0 ;
468 a10 = -dt*(R1*coeff(4) + R1*R1*coeff(5) + coeff(6))/R2 ;
469 a11 = -dt*(coeff(4) + 2*R1*coeff(5)) ;
470 a12 = -dt*R2*coeff(5) ;
471 a13 = 0 ;
472 a20 = -dt*R1/(R2*R2)*(coeff(1) + R1*coeff(2) + R1*R1*coeff(3)) ;
473 a21 = -dt/R2*(coeff(1) + 2*R1*coeff(2) + 3*R1*R1*coeff(3)) ;
474 a22 = -dt*(coeff(2) + 3*R1*coeff(3)) ;
475 a23 = -dt*R2*coeff(3) ;
476 a24 = 0 ;
477 type_dal = ((0.1*(fabs(a20)+fabs(a21)+fabs(a22)+fabs(a23))*nr*nr*nr
478 < 1.) ? O2DEGE_SMALL : O2DEGE_LARGE ) ;
479 break ;
480 case false:
481 a00 = R1*R1 - dt*(coeff(7)*R1*R1 + coeff(8)*R1 + coeff(9)) ;
482 a01 = 2*R1*R2 - dt*(2*R1*R2*coeff(7) + R2*coeff(8)) ;
483 a02 = R2*R2*(1 - dt*coeff(7)) ;
484 a10 = -dt*R1/R2*(R1*coeff(4) + R1*R1*coeff(5) + coeff(6)) ;
485 a11 = -dt*(2*R1*coeff(4) + 3*R1*R1*coeff(5) + coeff(6)) ;
486 a12 = -dt*(R2*coeff(4) + 3*R1*R2*coeff(5)) ;
487 a13 = -dt*R2*R2*coeff(5) ;
488 a20 = -dt*(R1*R1)/(R2*R2)*(coeff(1) + R1*coeff(2) + R1*R1*coeff(3)) ;
489 a21 = -dt*R1/R2*(2*coeff(1) + 3*R1*coeff(2) + 4*R1*R1*coeff(3)) ;
490 a22 = -dt*(coeff(1) + 3*R1*coeff(2) + 6*R1*R1*coeff(3)) ;
491 a23 = -dt*(R2*coeff(2) + 4*R1*R2*coeff(3)) ;
492 a24 = -dt*R2*R2*coeff(3) ;
493 type_dal = ((0.1*(fabs(a20)+fabs(a21)+fabs(a22)+fabs(a23)+fabs(a24))
494 *nr*nr*nr < 1.) ? O2NOND_SMALL : O2NOND_LARGE ) ;
495 break ;
496 }
497 if (fabs(a00)<1.e-15) a00 = 0 ;
498 if (fabs(a01)<1.e-15) a01 = 0 ;
499 if (fabs(a02)<1.e-15) a02 = 0 ;
500 if (fabs(a10)<1.e-15) a10 = 0 ;
501 if (fabs(a11)<1.e-15) a11 = 0 ;
502 if (fabs(a12)<1.e-15) a12 = 0 ;
503 if (fabs(a13)<1.e-15) a13 = 0 ;
504 if (fabs(a20)<1.e-15) a20 = 0 ;
505 if (fabs(a21)<1.e-15) a21 = 0 ;
506 if (fabs(a22)<1.e-15) a22 = 0 ;
507 if (fabs(a23)<1.e-15) a23 = 0 ;
508 if (fabs(a24)<1.e-15) a24 = 0 ;
509
510
511
512 Diff_id id(R_JACO02, nr) ;
513 if (fabs(a00)>1.e-15) {
514 operateur = a00*id ;
515 }
516 else{
517 operateur.set_etat_qcq() ;
518 for (int i=0; i<nr; i++)
519 for (int j=0; j<nr; j++)
520 operateur.set(i,j) = 0. ;
521 }
522 Diff_mx op01(R_JACO02, nr) ; const Matrice& m01 = op01.get_matrice() ;
523 Diff_mx2 op02(R_JACO02, nr) ; const Matrice& m02 = op02.get_matrice() ;
524 Diff_dsdx op10(R_JACO02, nr) ; const Matrice& m10 = op10.get_matrice() ;
525 Diff_xdsdx op11(R_JACO02, nr) ; const Matrice& m11 = op11.get_matrice() ;
526 Diff_x2dsdx op12(R_JACO02, nr) ; const Matrice& m12 = op12.get_matrice() ;
527 Diff_x3dsdx op13(R_JACO02, nr) ; const Matrice& m13 = op13.get_matrice() ;
528 Diff_dsdx2 op20(R_JACO02, nr) ; const Matrice& m20 = op20.get_matrice() ;
529 Diff_xdsdx2 op21(R_JACO02, nr) ; const Matrice& m21 = op21.get_matrice() ;
530 Diff_x2dsdx2 op22(R_JACO02, nr) ; const Matrice& m22 = op22.get_matrice() ;
531 Diff_x3dsdx2 op23(R_JACO02, nr) ; const Matrice& m23 = op23.get_matrice() ;
532 Diff_x4dsdx2 op24(R_JACO02, nr) ; const Matrice& m24 = op24.get_matrice() ;
533
534 for (int i=0; i<nr; i++) {
535 int jmin = (i>3 ? i-3 : 0) ;
536 int jmax = (i<nr-9 ? i+10 : nr) ;
537 for (int j=jmin ; j<jmax; j++)
538 operateur.set(i,j) += a01*m01(i,j) + a02*m02(i,j)
539 + a10*m10(i,j) + a11*m11(i,j) + a12*m12(i,j)
540 + a13*m13(i,j) + a20*m20(i,j) + a21*m21(i,j)
541 + a22*m22(i,j) + a23*m23(i,j) + a24*m24(i,j) ;
542
543 }
544}
545
546
547
548
549 //--------------------------
550 //- La routine a appeler ---
551 //----------------------------
552void get_operateur_dal(const Param& par, const int& lzone,
553 const int& base_r, int& type_dal, Matrice& operateur)
554{
555
556 // Routines de derivation
557 static void (*get_operateur_dal[MAX_BASE])(const Param&,
558 const int&, int&, Matrice&) ;
559 static int nap = 0 ;
560
561 // Premier appel
562 if (nap==0) {
563 nap = 1 ;
564 for (int i=0 ; i<MAX_BASE ; i++)
565 get_operateur_dal[i] = _get_operateur_dal_pas_prevu ;
566
567 // Les routines existantes
568 get_operateur_dal[R_CHEB >> TRA_R] = _get_operateur_dal_r_cheb ;
569 get_operateur_dal[R_CHEBP >> TRA_R] = _get_operateur_dal_r_chebp ;
570 get_operateur_dal[R_CHEBI >> TRA_R] = _get_operateur_dal_r_chebi ;
571 get_operateur_dal[R_JACO02 >> TRA_R] = _get_operateur_dal_r_jaco02 ;
572 }
573
574 get_operateur_dal[base_r](par, lzone, type_dal, operateur) ;
575}
576}
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:263
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:67