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