LORENE
ope_helmholtz_minus_2d_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_2d_mat.C,v 1.4 2016/12/05 16:18:11 j_novak Exp $
25  * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_helmholtz_minus_2d/ope_helmholtz_minus_2d_mat.C,v 1.4 2016/12/05 16:18:11 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_2d_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_2d_mat_r_chebu_deux(int,int,double, double) ;
55 
56 Matrice _helmholtz_minus_2d_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_2d_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_2d_mat_r_chebu_deux (int n, int l, double masse,
73  double alpha) {
74 
75  Matrice res(n-2, n-2) ;
76  res.set_etat_qcq() ;
77  double* vect = new double[n] ;
78  double* vect_bis = new double[n] ;
79  double* vect_dd = new double[n] ;
80  double* vect_d = new double[n] ;
81 
82  for (int i=0 ; i<n-2 ; i++) {
83  for (int j=0 ; j<n ; j++)
84  vect[j] = 0 ;
85  vect[i] = 2*i+3 ;
86  vect[i+1] = -4*i-4 ;
87  vect[i+2] = 2*i+1 ;
88 
89  // Der sec.
90  for (int j=0 ; j<n ; j++)
91  vect_bis[j] = vect[j] ;
92 
93  d2sdx2_1d (n, &vect_bis, R_CHEBU) ; // appel dans le cas unsurr
94  mult2_xm1_1d_cheb (n, vect_bis, vect_dd) ; // multiplication par (x-1)^2
95 
96  // Der simple
97  for (int j=0 ; j<n ; j++)
98  vect_bis[j] = vect[j] ;
99 
100  dsdx_1d (n, &vect_bis, R_CHEBU) ; // appel dans le cas unsurr
101  mult_xm1_1d_cheb (n, vect_bis, vect_d) ; // multiplication par (x-1)
102 
103  // Mass term
104  for (int j=0 ; j<n ; j++)
105  vect_bis[j] = vect[j] ;
106  sx2_1d (n, &vect_bis, R_CHEBU) ;
107 
108  for (int j=0 ; j<n-2 ; j++)
109  res.set(j,i) = vect_dd[j] + vect_d[j] - l*l*vect[j] - masse*masse/alpha/alpha*vect_bis[j] ;
110  }
111 
112  delete [] vect ;
113  delete [] vect_bis ;
114  delete [] vect_dd ;
115  delete [] vect_d ;
116 
117  return res ;
118 }
119 
120 
121  //-------------------------
122  //-- CAS R_CHEB -----
123  //-----------------------
124 
125 
126 Matrice _helmholtz_minus_2d_mat_r_cheb (int n, int l, double masse,
127  double alf, double bet, int) {
128 
129  double echelle = bet / alf ;
130 
131  Matrice dd(n, n) ;
132  dd.set_etat_qcq() ;
133  Matrice xd(n, n) ;
134  xd.set_etat_qcq() ;
135  Matrice xx(n, n) ;
136  xx.set_etat_qcq() ;
137 
138  double* vect = new double[n] ;
139 
140  for (int i=0 ; i<n ; i++) {
141  for (int j=0 ; j<n ; j++)
142  vect[j] = 0 ;
143  vect[i] = 1 ;
144  d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
145  vect[i] -= masse*masse*alf*alf ;
146  for (int j=0 ; j<n ; j++)
147  dd.set(j, i) = vect[j]*echelle*echelle ;
148  }
149 
150  for (int i=0 ; i<n ; i++) {
151  for (int j=0 ; j<n ; j++)
152  vect[j] = 0 ;
153  vect[i] = 1 ;
154  d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
155  vect[i] -= masse*masse*alf*alf ;
156  multx_1d (n, &vect, R_CHEB) ;
157  for (int j=0 ; j< n ; j++)
158  dd.set(j, i) += 2*echelle*vect[j] ;
159  }
160 
161  for (int i=0 ; i<n ; i++) {
162  for (int j=0 ; j<n ; j++)
163  vect[j] = 0 ;
164  vect[i] = 1 ;
165  d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
166  vect[i] -= masse*masse*alf*alf ;
167  multx_1d (n, &vect, R_CHEB) ;
168  multx_1d (n, &vect, R_CHEB) ;
169  for (int j=0 ; j<n ; j++)
170  dd.set(j, i) += vect[j] ;
171  }
172 
173  for (int i=0 ; i<n ; i++) {
174  for (int j=0 ; j<n ; j++)
175  vect[j] = 0 ;
176  vect[i] = 1 ;
177  sxdsdx_1d (n, &vect, R_CHEB) ;
178  for (int j=0 ; j<n ; j++)
179  xd.set(j, i) = vect[j]*echelle ;
180  }
181 
182  for (int i=0 ; i<n ; i++) {
183  for (int j=0 ; j<n ; j++)
184  vect[j] = 0 ;
185  vect[i] = 1 ;
186  sxdsdx_1d (n, &vect, R_CHEB) ;
187  multx_1d (n, &vect, R_CHEB) ;
188  for (int j=0 ; j<(n>i+1 ? i+1 : n) ; j++)
189  xd.set(j, i) += vect[j] ;
190  }
191 
192  for (int i=0 ; i<n ; i++) {
193  for (int j=0 ; j<n ; j++)
194  vect[j] = 0 ;
195  vect[i] = 1 ;
196  sx2_1d (n, &vect, R_CHEB) ;
197  for (int j=0 ; j<n ; j++)
198  xx.set(j, i) = vect[j] ;
199  }
200 
201  delete [] vect ;
202 
203  Matrice res(n, n) ;
204  res = dd+xd-l*l*xx ;
205 
206  return res ;
207 }
208 
209 
211  if (ope_mat != 0x0)
212  delete ope_mat ;
213 
214  // Routines de derivation
215  static Matrice (*helmholtz_minus_2d_mat[MAX_BASE])(int, int, double,
216  double, double, int);
217  static int nap = 0 ;
218 
219  // Premier appel
220  if (nap==0) {
221  nap = 1 ;
222  for (int i=0 ; i<MAX_BASE ; i++) {
223  helmholtz_minus_2d_mat[i] = _helmholtz_minus_2d_mat_pas_prevu ;
224  }
225  // Les routines existantes
226  helmholtz_minus_2d_mat[R_CHEB >> TRA_R] = _helmholtz_minus_2d_mat_r_cheb ;
227  helmholtz_minus_2d_mat[R_CHEBU >> TRA_R] = _helmholtz_minus_2d_mat_r_chebu ;
228  }
229  ope_mat = new Matrice(helmholtz_minus_2d_mat[base_r](nr, l_quant, masse,
230  alpha, beta, dzpuis)) ;
231 }
232 }
double alpha
Parameter of the associated mapping.
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.
double masse
The mass term.
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
Matrix handling.
Definition: matrice.h:152
int dzpuis
the associated dzpuis, if in the compactified domain.
int nr
Number of radial points.
virtual void do_ope_mat() const
Computes the matrix of the operator.
#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
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166