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 
108 namespace Lorene {
109 void _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 
117 void _get_operateur_dal_r_cheb(const Param& par, const int& lz,
118 int& 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 
234 void _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 
331 void _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 
429 void _get_operateur_dal_r_jaco02(const Param& par, const int& lz,
430 int& 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  //----------------------------
552 void 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 }
O2DEGE_LARGE
#define O2DEGE_LARGE
Operateur du deuxieme ordre degenere .
Definition: type_parite.h:282
Lorene
Lorene prototypes.
Definition: app_hor.h:67
O2NOND_LARGE
#define O2NOND_LARGE
Operateur du deuxieme ordre non degenere .
Definition: type_parite.h:286
R_JACO02
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
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_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
ORDRE1_SMALL
#define ORDRE1_SMALL
Operateur du premier ordre, .
Definition: type_parite.h:276
O2DEGE_SMALL
#define O2DEGE_SMALL
Operateur du deuxieme ordre degenere .
Definition: type_parite.h:280
ORDRE1_LARGE
#define ORDRE1_LARGE
Operateur du premier ordre .
Definition: type_parite.h:278
MAX_BASE
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166