LORENE
ope_sec_order_mat.C
1 /*
2  * Copyright (c) 2003 Philippe Grandclement
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 version 2
8  * as published by the Free Software Foundation.
9  *
10  * LORENE is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with LORENE; if not, write to the Free Software
17  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18  *
19  */
20 
21 
22 
23 /*
24  * $Id: ope_sec_order_mat.C,v 1.4 2016/12/05 16:18:13 j_novak Exp $
25  * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_sec_order/ope_sec_order_mat.C,v 1.4 2016/12/05 16:18:13 j_novak Exp $
26  *
27  */
28 #include <cmath>
29 #include <cstdlib>
30 
31 #include "proto.h"
32 #include "ope_elementary.h"
33 
34 
35  //-----------------------------------
36  // Routine pour les cas non prevus --
37  //-----------------------------------
38 
39 namespace Lorene {
40 Matrice _sec_order_mat_pas_prevu(int, double, double, double, double) {
41  cout << "Sec_order : base not implemented..." << endl ;
42  abort() ;
43  exit(-1) ;
44  Matrice res(1, 1) ;
45  return res;
46 }
47 
48 
49 
50  //-------------------------
51  //-- CAS R_CHEB -----
52  //------------------------
53 
54 Matrice _sec_order_mat_r_cheb (int n, double alpha,
55  double a, double b, double c) {
56 
57  double* vect = new double[n] ;
58 
59  Matrice dd (n,n) ;
60  dd.set_etat_qcq() ;
61  Matrice df (n,n) ;
62  df.set_etat_qcq() ;
63  Matrice ff (n,n) ;
64  ff.set_etat_qcq() ;
65 
66 
67  // Calcul terme en d^2...
68  for (int i=0 ; i<n ; i++) {
69  for (int j=0 ; j<n ; j++)
70  vect[j] = 0 ;
71  vect[i] = 1 ;
72 
73  d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
74 
75  for (int j=0 ; j<n ; j++)
76  dd.set(j,i) = vect[j] ;
77  }
78 
79  // Calcul terme en d...
80  for (int i=0 ; i<n ; i++) {
81  for (int j=0 ; j<n ; j++)
82  vect[j] = 0 ;
83  vect[i] = 1 ;
84 
85  sxdsdx_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
86 
87  for (int j=0 ; j<n ; j++)
88  df.set(j,i) = vect[j] ;
89  }
90 
91  // Calcul terme Id
92  for (int i=0 ; i<n ; i++) {
93  for (int j=0 ; j<n ; j++)
94  ff.set(j,i) = 0 ;
95  ff.set(i,i) = 1 ;
96  }
97 
98  delete [] vect ;
99 
100  return a*dd/alpha/alpha+b*df/alpha+c*ff ;
101 }
102 
104  if (ope_mat != 0x0)
105  delete ope_mat ;
106 
107  // Routines de derivation
108  static Matrice (*sec_order_mat[MAX_BASE])(int, double,
109  double, double, double);
110  static int nap = 0 ;
111 
112  // Premier appel
113  if (nap==0) {
114  nap = 1 ;
115  for (int i=0 ; i<MAX_BASE ; i++) {
116  sec_order_mat[i] = _sec_order_mat_pas_prevu ;
117  }
118  // Les routines existantes
119  sec_order_mat[R_CHEB >> TRA_R] = _sec_order_mat_r_cheb ;
120  }
121  ope_mat = new Matrice(sec_order_mat[base_r](nr, alpha,
122  a_param, b_param, c_param)) ;
123 }
124 }
double alpha
Parameter of the associated mapping.
Lorene prototypes.
Definition: app_hor.h:67
Matrice * ope_mat
Pointer on the matrix representation of the operator.
int base_r
Radial basis of decomposition.
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
double b_param
The parameter .
double c_param
The parameter .
Matrix handling.
Definition: matrice.h:152
double a_param
The parameter .
virtual void do_ope_mat() const
Computes the matrix of the operator.
int nr
Number of radial points.
#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