LORENE
helmholtz_minus_mat.C
1/*
2 * Copyright (c) 1999-2001 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 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
23char helmholtz_minus_mat_C[] = "$$" ;
24
25/*
26 * $Id: helmholtz_minus_mat.C,v 1.8 2014/10/13 08:53:28 j_novak Exp $
27 * $Log: helmholtz_minus_mat.C,v $
28 * Revision 1.8 2014/10/13 08:53:28 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.7 2014/10/06 15:16:08 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.6 2008/07/09 06:51:58 p_grandclement
35 * some corrections to helmholtz minus in the nucleus
36 *
37 * Revision 1.5 2008/07/08 11:45:28 p_grandclement
38 * Add helmholtz_minus in the nucleus
39 *
40 * Revision 1.4 2004/08/24 09:14:44 p_grandclement
41 * Addition of some new operators, like Poisson in 2d... It now requieres the
42 * GSL library to work.
43 *
44 * Also, the way a variable change is stored by a Param_elliptic is changed and
45 * no longer uses Change_var but rather 2 Scalars. The codes using that feature
46 * will requiere some modification. (It should concern only the ones about monopoles)
47 *
48 * Revision 1.3 2004/01/15 09:15:37 p_grandclement
49 * Modification and addition of the Helmholtz operators
50 *
51 * Revision 1.2 2003/12/11 15:57:26 p_grandclement
52 * include stdlib.h encore ...
53 *
54 * Revision 1.1 2003/12/11 14:48:49 p_grandclement
55 * Addition of ALL (and that is a lot !) the files needed for the general elliptic solver ... UNDER DEVELOPEMENT...
56 *
57 *
58 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/helmholtz_minus_mat.C,v 1.8 2014/10/13 08:53:28 j_novak Exp $
59 *
60 */
61#include <cstdlib>
62
63#include "matrice.h"
64#include "type_parite.h"
65#include "proto.h"
66#include "diff.h"
67
68 //-----------------------------------
69 // Routine pour les cas non prevus --
70 //-----------------------------------
71
72namespace Lorene {
73Matrice _helmholtz_minus_mat_pas_prevu(int, int, double, double, double) {
74 cout << "Helmholtz minus : base not implemented..." << endl ;
75 abort() ;
76 exit(-1) ;
77 Matrice res(1, 1) ;
78 return res;
79}
80
81 //-------------------------
82 //-- CAS R_CHEBU -----
83 //------------------------
84
85Matrice _helmholtz_minus_mat_r_chebu (int n, int lq, double alpha,
86 double, double masse) {
87
88 assert (masse > 0) ;
89
90 Matrice res(n-2, n-2) ;
91 res.set_etat_qcq() ;
92
93 double* vect = new double[n] ;
94 double* vect_bis = new double[n] ;
95 double* vect_dd = new double[n] ;
96
97 for (int i=0 ; i<n-2 ; i++) {
98
99 for (int j=0 ; j<n ; j++)
100 vect[j] = 0 ;
101 vect[i] = 2*i+3 ;
102 vect[i+1] = -4*i-4 ;
103 vect[i+2] = 2*i+1 ;
104
105 for (int j=0 ; j<n ; j++)
106 vect_bis[j] = vect[j] ;
107
108 d2sdx2_1d (n, &vect_bis, R_CHEBU) ; // appel dans le cas unsurr
109 mult2_xm1_1d_cheb (n, vect_bis, vect_dd) ; // multiplication par (x-1)^2
110
111 // Mass term
112 for (int j=0 ; j<n ; j++)
113 vect_bis[j] = vect[j] ;
114 sx2_1d (n, &vect_bis, R_CHEBU) ;
115
116 for (int j=0 ; j<n-2 ; j++)
117 res.set(j,i) = vect_dd[j] - lq*(lq+1)*vect[j]
118 - masse*masse*vect_bis[j]/alpha/alpha ;
119 }
120
121 delete [] vect ;
122 delete [] vect_bis ;
123 delete [] vect_dd ;
124
125 return res ;
126}
127
128
129 //-------------------------
130 //-- CAS R_CHEB -----
131 //------------------------
132
133Matrice _helmholtz_minus_mat_r_cheb (int n, int lq, double alpha, double beta,
134 double masse) {
135
136 assert (masse > 0) ;
137
138 double echelle = beta / alpha ;
139
140 Matrice dd(n, n) ;
141 dd.set_etat_qcq() ;
142 Matrice xd(n, n) ;
143 xd.set_etat_qcq() ;
144 Matrice xx(n, n) ;
145 xx.set_etat_qcq() ;
146
147 double* vect = new double[n] ;
148
149 for (int i=0 ; i<n ; i++) {
150 for (int j=0 ; j<n ; j++)
151 vect[j] = 0 ;
152 vect[i] = 1 ;
153 d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
154 vect[i] -= masse*masse*alpha*alpha ;
155 for (int j=0 ; j<n ; j++)
156 dd.set(j, i) = vect[j]*echelle*echelle ;
157 }
158
159 for (int i=0 ; i<n ; i++) {
160 for (int j=0 ; j<n ; j++)
161 vect[j] = 0 ;
162 vect[i] = 1 ;
163 d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
164 vect[i] -= masse*masse*alpha*alpha ;
165 multx_1d (n, &vect, R_CHEB) ;
166 for (int j=0 ; j<n ; j++)
167 dd.set(j, i) += 2*echelle*vect[j] ;
168 }
169
170 for (int i=0 ; i<n ; i++) {
171 for (int j=0 ; j<n ; j++)
172 vect[j] = 0 ;
173 vect[i] = 1 ;
174 d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
175 vect[i] -= masse*masse*alpha*alpha ;
176 multx_1d (n, &vect, R_CHEB) ;
177 multx_1d (n, &vect, R_CHEB) ;
178 for (int j=0 ; j<n ; j++)
179 dd.set(j, i) += vect[j] ;
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 for (int j=0 ; j<n ; j++)
188 xd.set(j, i) = vect[j]*echelle ;
189 }
190
191 for (int i=0 ; i<n ; i++) {
192 for (int j=0 ; j<n ; j++)
193 vect[j] = 0 ;
194 vect[i] = 1 ;
195 sxdsdx_1d (n, &vect, R_CHEB) ;
196 multx_1d (n, &vect, R_CHEB) ;
197 for (int j=0 ; j<n ; j++)
198 xd.set(j, i) += vect[j] ;
199 }
200
201 for (int i=0 ; i<n ; i++) {
202 for (int j=0 ; j<n ; j++)
203 vect[j] = 0 ;
204 vect[i] = 1 ;
205 sx2_1d (n, &vect, R_CHEB) ;
206 for (int j=0 ; j<n ; j++)
207 xx.set(j, i) = vect[j] ;
208 }
209
210 delete [] vect ;
211
212 Matrice res(n, n) ;
213 res = dd+2*xd - lq*(lq+1)*xx;
214
215 return res ;
216}
217
218
219 //-------------------------
220 //-- CAS R_CHEBP -----
221 //--------------------------
222
223
224Matrice _helmholtz_minus_mat_r_chebp (int n, int lq, double alpha, double, double masse) {
225
226 if (lq==0) {
227 Diff_dsdx2 d2(R_CHEBP, n) ;
228 Diff_sxdsdx sxd(R_CHEBP, n) ;
229 Diff_id xx (R_CHEBP, n) ;
230
231 return Matrice(d2 + 2.*sxd -masse*masse*alpha*alpha*xx) ;
232 }
233 else {
234 Matrice res(n-1, n-1) ;
235 res.set_etat_qcq() ;
236
237 double* vect = new double[n] ;
238
239 double* vect_sx2 = new double[n] ;
240 double* vect_sxd = new double[n] ;
241 double* vect_dd = new double[n] ;
242
243 for (int i=0 ; i<n-1 ; i++) {
244 for (int j=0 ; j<n ; j++)
245 vect[j] = 0 ;
246 vect[i] = 1. ;
247 vect[i+1] = 1. ;
248
249
250 for (int j=0 ; j<n ; j++)
251 vect_dd[j] = vect[j] ;
252 d2sdx2_1d (n, &vect_dd, R_CHEBP) ; // appel dans le cas chebp
253 for (int j=0 ; j<n ; j++)
254 vect_sxd[j] = vect[j] ;
255 sxdsdx_1d (n, &vect_sxd, R_CHEBP) ; // appel dans le cas chebp
256 for (int j=0 ; j<n ; j++)
257 vect_sx2[j] = vect[j] ;
258 sx2_1d (n, &vect_sx2, R_CHEBP) ; // appel dans le cas chebp
259
260 for (int j=0 ; j<n-1 ; j++)
261 res.set(j,i) = vect_dd[j] +2*vect_sxd[j] - lq*(lq+1)*vect_sx2[j] - masse*masse*alpha*alpha*vect[j] ;
262 }
263 delete [] vect ;
264 delete [] vect_sx2 ;
265 delete [] vect_sxd ;
266 delete [] vect_dd ;
267
268 return res ;
269 }
270}
271
272
273
274 //------------------------
275 //-- CAS R_CHEBI ----
276 //------------------------
277
278
279Matrice _helmholtz_minus_mat_r_chebi (int n, int lq, double alpha, double, double masse) {
280
281 if (lq==1) {
282 Diff_dsdx2 d2(R_CHEBI, n) ;
283 Diff_sxdsdx sxd(R_CHEBI, n) ;
284 Diff_sx2 sx2(R_CHEBI, n) ;
285 Diff_id xx(R_CHEBI, n) ;
286
287 return Matrice(d2 + 2.*sxd - (lq*(lq+1))*sx2- masse*masse*alpha*alpha*xx) ;
288 }
289 else {
290 Matrice res(n-1, n-1) ;
291 res.set_etat_qcq() ;
292
293 double* vect = new double[n] ;
294
295 double* vect_sx2 = new double[n] ;
296 double* vect_sxd = new double[n] ;
297 double* vect_dd = new double[n] ;
298
299 for (int i=0 ; i<n-1 ; i++) {
300 for (int j=0 ; j<n ; j++)
301 vect[j] = 0 ;
302 vect[i] = (2*i+3) ;
303 vect[i+1] = (2*i+1) ;
304
305
306 for (int j=0 ; j<n ; j++)
307 vect_dd[j] = vect[j] ;
308 d2sdx2_1d (n, &vect_dd, R_CHEBI) ; // appel dans le cas chebi
309 for (int j=0 ; j<n ; j++)
310 vect_sxd[j] = vect[j] ;
311 sxdsdx_1d (n, &vect_sxd, R_CHEBI) ; // appel dans le cas chebi
312 for (int j=0 ; j<n ; j++)
313 vect_sx2[j] = vect[j] ;
314 sx2_1d (n, &vect_sx2, R_CHEBI) ; // appel dans le cas chebi
315
316 for (int j=0 ; j<n-1 ; j++)
317 res.set(j,i) = vect_dd[j] +2*vect_sxd[j] - lq*(lq+1)*vect_sx2[j] - masse*masse*alpha*alpha*vect[j] ;
318 }
319 delete [] vect ;
320 delete [] vect_sx2 ;
321 delete [] vect_sxd ;
322 delete [] vect_dd ;
323
324 return res ;
325 }
326}
327
328
329
330 //--------------------------
331 //- La routine a appeler ---
332 //----------------------------
333
334Matrice helmholtz_minus_mat(int n, int lq,
335 double alpha, double beta, double masse,
336 int base_r)
337{
338
339 // Routines de derivation
340 static Matrice (*helmholtz_minus_mat[MAX_BASE])(int, int,
341 double, double, double);
342 static int nap = 0 ;
343
344 // Premier appel
345 if (nap==0) {
346 nap = 1 ;
347 for (int i=0 ; i<MAX_BASE ; i++) {
348 helmholtz_minus_mat[i] = _helmholtz_minus_mat_pas_prevu ;
349 }
350 // Les routines existantes
351 helmholtz_minus_mat[R_CHEB >> TRA_R] = _helmholtz_minus_mat_r_cheb ;
352 helmholtz_minus_mat[R_CHEBU >> TRA_R] = _helmholtz_minus_mat_r_chebu ;
353 helmholtz_minus_mat[R_CHEBP >> TRA_R] = _helmholtz_minus_mat_r_chebp ;
354 helmholtz_minus_mat[R_CHEBI >> TRA_R] = _helmholtz_minus_mat_r_chebi ;
355 }
356
357 Matrice res(helmholtz_minus_mat[base_r](n, lq, alpha, beta, masse)) ;
358 return res ;
359}
360
361}
Lorene::Diff_dsdx2
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:172
Lorene::Diff_id
Class for the elementary differential operator Identity (see the base class Diff ).
Definition diff.h:210
Lorene::Diff_sx2
Class for the elementary differential operator division by (see the base class Diff ).
Definition diff.h:369
Lorene::Diff_sxdsdx
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:450
Lorene::Matrice
Matrix handling.
Definition matrice.h:152
MAX_BASE
#define MAX_BASE
Nombre max. de bases differentes.
Definition type_parite.h:144
R_CHEBU
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition type_parite.h:180
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition type_parite.h:170
TRA_R
#define TRA_R
Translation en R, used for a bitwise shift (in hex).
Definition type_parite.h:158
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition type_parite.h:166
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
Lorene
Lorene prototypes.
Definition app_hor.h:64