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
21char ope_helmholtz_minus_2d_mat_C[] = "$Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_helmholtz_minus_2d/ope_helmholtz_minus_2d_mat.C,v 1.3 2014/10/13 08:53:33 j_novak Exp $" ;
22
23/*
24 * $Id: ope_helmholtz_minus_2d_mat.C,v 1.3 2014/10/13 08:53:33 j_novak Exp $
25 * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_helmholtz_minus_2d/ope_helmholtz_minus_2d_mat.C,v 1.3 2014/10/13 08:53:33 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
38namespace Lorene {
39Matrice _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
54Matrice _helmholtz_minus_2d_mat_r_chebu_deux(int,int,double, double) ;
55
56Matrice _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
72Matrice _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
126Matrice _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
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 }
230 alpha, beta, dzpuis)) ;
231}
232}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Matrix handling.
Definition matrice.h:152
Matrice * ope_mat
Pointer on the matrix representation of the operator.
double beta
Parameter of the associated mapping.
double alpha
Parameter of the associated mapping.
int base_r
Radial basis of decomposition.
int nr
Number of radial points.
virtual void do_ope_mat() const
Computes the matrix of the operator.
int dzpuis
the associated dzpuis, if in the compactified domain.
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
#define R_CHEB
base de Chebychev ordinaire (fin)
Lorene prototypes.
Definition app_hor.h:64