LORENE
ope_helmholtz_minus_pseudo_1d_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_helmholtz_minus_pseudo_1d_mat.C,v 1.4 2016/12/05 16:18:12 j_novak Exp $
25  * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_helmholtz_minus_pseudo_1d/ope_helmholtz_minus_pseudo_1d_mat.C,v 1.4 2016/12/05 16:18:12 j_novak Exp $
26  *
27  */
28 #include <cmath>
29 #include <cstdlib>
30 
31 #include "proto.h"
32 #include "ope_elementary.h"
33 
34  //-----------------------------------
35  // Routine pour les cas non prevus --
36  //-----------------------------------
37 
38 namespace Lorene {
39 Matrice _helmholtz_minus_pseudo_1d_mat_pas_prevu(int, int, double, double,
40  double, int) {
41  cout << "Operateur pas prevu..." << endl ;
42  abort() ;
43  exit(-1) ;
44  Matrice res(1, 1) ;
45  return res;
46 }
47 
48 
49 
50  //-------------------------
51  //-- CAS R_CHEBU -----
52  //-------------------------
53 
54 Matrice _helmholtz_minus_pseudo_1d_mat_r_chebu_deux(int,int,double, double) ;
55 
56 Matrice _helmholtz_minus_pseudo_1d_mat_r_chebu( int n, int l, double masse,
57  double alpha, double, int puis) {
58  Matrice res(n-2, n-2) ;
59  res.set_etat_qcq() ;
60  switch (puis) {
61  case 2 :
62  res = _helmholtz_minus_pseudo_1d_mat_r_chebu_deux (n, l,masse, alpha) ;
63  break ;
64  default :
65  abort() ;
66  exit(-1) ;
67  }
68  return res ;
69 }
70 
71  //Cas ou dzpuis = 2
72 Matrice _helmholtz_minus_pseudo_1d_mat_r_chebu_deux (int n, int l, double masse,
73  double alpha) {
74 
75 
76  Matrice res(n-2, n-2) ;
77  res.set_etat_qcq() ;
78  double* vect = new double[n] ;
79  double* vect_bis = new double[n] ;
80  double* vect_dd = new double[n] ;
81  double* vect_d = new double[n] ;
82 
83  for (int i=0 ; i<n-2 ; i++) {
84  for (int j=0 ; j<n ; j++)
85  vect[j] = 0 ;
86  vect[i] = 2*i+3 ;
87  vect[i+1] = -4*i-4 ;
88  vect[i+2] = 2*i+1 ;
89 
90  // Der sec.
91  for (int j=0 ; j<n ; j++)
92  vect_bis[j] = vect[j] ;
93 
94  d2sdx2_1d (n, &vect_bis, R_CHEBU) ; // appel dans le cas unsurr
95  mult2_xm1_1d_cheb (n, vect_bis, vect_dd) ; // multiplication par (x-1)^2
96 
97  // Der simple
98  for (int j=0 ; j<n ; j++)
99  vect_bis[j] = vect[j] ;
100 
101  dsdx_1d (n, &vect_bis, R_CHEBU) ; // appel dans le cas unsurr
102  mult_xm1_1d_cheb (n, vect_bis, vect_d) ; // multiplication par (x-1)
103 
104  // Mass term
105  for (int j=0 ; j<n ; j++)
106  vect_bis[j] = vect[j] ;
107  sx2_1d (n, &vect_bis, R_CHEBU) ;
108 
109  for (int j=0 ; j<n-2 ; j++)
110  res.set(j,i) = vect_dd[j] + 2*vect_d[j] - l*(l-1)*vect[j] - masse*masse/alpha/alpha*vect_bis[j] ;
111  }
112 
113  delete [] vect ;
114  delete [] vect_bis ;
115  delete [] vect_dd ;
116 
117  return res ;
118 }
119 
120 
121 
123  if (ope_mat != 0x0)
124  delete ope_mat ;
125 
126  // Routines de derivation
127  static Matrice (*helmholtz_minus_pseudo_1d_mat[MAX_BASE])(int, int, double,
128  double, double, int);
129  static int nap = 0 ;
130 
131  // Premier appel
132  if (nap==0) {
133  nap = 1 ;
134  for (int i=0 ; i<MAX_BASE ; i++) {
135  helmholtz_minus_pseudo_1d_mat[i] = _helmholtz_minus_pseudo_1d_mat_pas_prevu ;
136  }
137  // Les routines existantes
138  helmholtz_minus_pseudo_1d_mat[R_CHEBU >> TRA_R] = _helmholtz_minus_pseudo_1d_mat_r_chebu ;
139  }
140  ope_mat = new Matrice(helmholtz_minus_pseudo_1d_mat[base_r](nr, l_quant, masse,
141  alpha, beta, dzpuis)) ;
142 }
143 }
double alpha
Parameter of the associated mapping.
int dzpuis
the associated dzpuis, if in the compactified domain.
double beta
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
Matrix handling.
Definition: matrice.h:152
virtual void do_ope_mat() const
Computes the matrix of the operator.
int nr
Number of radial points.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144