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 char 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 
105 namespace Lorene {
106 void _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 
114 void _get_operateur_dal_r_cheb(const Param& par, const int& lz,
115 int& 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 
231 void _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 
328 void _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 
426 void _get_operateur_dal_r_jaco02(const Param& par, const int& lz,
427 int& 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  //----------------------------
549 void 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 }
#define O2DEGE_LARGE
Operateur du deuxieme ordre degenere .
Definition: type_parite.h:282
Lorene prototypes.
Definition: app_hor.h:64
#define O2NOND_LARGE
Operateur du deuxieme ordre non degenere .
Definition: type_parite.h:286
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
#define O2NOND_SMALL
Operateur du deuxieme ordre non degenere .
Definition: type_parite.h:284
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
#define ORDRE1_SMALL
Operateur du premier ordre, .
Definition: type_parite.h:276
#define O2DEGE_SMALL
Operateur du deuxieme ordre degenere .
Definition: type_parite.h:280
#define ORDRE1_LARGE
Operateur du premier ordre .
Definition: type_parite.h:278
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166