LORENE
helmholtz_plus_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_plus_mat_C[] = "$$" ;
24
25/*
26 * $Id: helmholtz_plus_mat.C,v 1.6 2014/10/13 08:53:28 j_novak Exp $
27 * $Log: helmholtz_plus_mat.C,v $
28 * Revision 1.6 2014/10/13 08:53:28 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.5 2014/10/06 15:16:08 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.4 2007/05/06 10:48:12 p_grandclement
35 * Modification of a few operators for the vorton project
36 *
37 * Revision 1.3 2004/01/15 09:15:37 p_grandclement
38 * Modification and addition of the Helmholtz operators
39 *
40 * Revision 1.2 2003/12/11 15:57:26 p_grandclement
41 * include stdlib.h encore ...
42 *
43 * Revision 1.1 2003/12/11 14:48:49 p_grandclement
44 * Addition of ALL (and that is a lot !) the files needed for the general elliptic solver ... UNDER DEVELOPEMENT...
45 *
46 *
47 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/helmholtz_plus_mat.C,v 1.6 2014/10/13 08:53:28 j_novak Exp $
48 *
49 */
50#include <cstdlib>
51
52#include "matrice.h"
53#include "type_parite.h"
54#include "proto.h"
55
56 //-----------------------------------
57 // Routine pour les cas non prevus --
58 //-----------------------------------
59
60namespace Lorene {
61Matrice _helmholtz_plus_mat_pas_prevu(int, int, double, double, double) {
62 cout << "Helmholtz plus : base not implemented..." << endl ;
63 abort() ;
64 exit(-1) ;
65 Matrice res(1, 1) ;
66 return res;
67}
68
69
70
71 //-------------------------
72 //-- CAS R_CHEBP -----
73 //--------------------------
74
75
76Matrice _helmholtz_plus_mat_r_chebp (int n, int lq, double alpha, double,
77 double masse) {
78 assert (masse > 0) ;
79
80 Matrice dd(n, n) ;
81 dd.set_etat_qcq() ;
82 Matrice xd(n, n) ;
83 xd.set_etat_qcq() ;
84 Matrice xx(n, n) ;
85 xx.set_etat_qcq() ;
86 Matrice sx2(n, n) ;
87 sx2.set_etat_qcq() ;
88
89 double* vect = new double[n] ;
90
91 for (int i=0 ; i<n ; i++) {
92 for (int j=0 ; j<n ; j++)
93 vect[j] = 0 ;
94 vect[i] = 1 ;
95 d2sdx2_1d (n, &vect, R_CHEBP) ;
96 for (int j=0 ; j<n ; j++)
97 dd.set(j, i) = vect[j] ;
98 }
99
100 for (int i=0 ; i<n ; i++) {
101 for (int j=0 ; j<n ; j++)
102 vect[j] = 0 ;
103 vect[i] = 1 ;
104 sxdsdx_1d (n, &vect, R_CHEBP) ;
105 for (int j=0 ; j<n ; j++)
106 xd.set(j, i) = vect[j] ;
107 }
108
109 for (int i=0 ; i<n ; i++) {
110 for (int j=0 ; j<n ; j++)
111 vect[j] = 0 ;
112 vect[i] = 1 ;
113 sx2_1d (n, &vect, R_CHEBP) ;
114 for (int j=0 ; j<n ; j++)
115 sx2.set(j, i) = vect[j] ;
116 }
117
118 for (int i=0 ; i<n ; i++) {
119 for (int j=0 ; j<n ; j++)
120 xx.set(i,j) = 0 ;
121 xx.set(i, i) = 1 ;
122 }
123
124 delete [] vect ;
125
126 Matrice res(n, n) ;
127 res = dd+2*xd-lq*(lq+1)*sx2+masse*masse*alpha*alpha*xx ;
128
129 return res ;
130}
131
132
133 //-------------------------
134 //-- CAS R_CHEB -----
135 //------------------------
136
137Matrice _helmholtz_plus_mat_r_cheb (int n, int lq, double alpha, double beta,
138 double masse) {
139
140
141
142 assert (masse > 0) ;
143
144 double echelle = beta / alpha ;
145
146 Matrice dd(n, n) ;
147 dd.set_etat_qcq() ;
148 Matrice xd(n, n) ;
149 xd.set_etat_qcq() ;
150 Matrice xx(n, n) ;
151 xx.set_etat_qcq() ;
152
153 double* vect = new double[n] ;
154
155 for (int i=0 ; i<n ; i++) {
156 for (int j=0 ; j<n ; j++)
157 vect[j] = 0 ;
158 vect[i] = 1 ;
159 d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
160 vect[i] += masse*masse*alpha*alpha ;
161 for (int j=0 ; j<n ; j++)
162 dd.set(j, i) = vect[j]*echelle*echelle ;
163 }
164
165 for (int i=0 ; i<n ; i++) {
166 for (int j=0 ; j<n ; j++)
167 vect[j] = 0 ;
168 vect[i] = 1 ;
169 d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
170 vect[i] += masse*masse*alpha*alpha ;
171 multx_1d (n, &vect, R_CHEB) ;
172 for (int j=0 ; j<n ; j++)
173 dd.set(j, i) += 2*echelle*vect[j] ;
174 }
175
176 for (int i=0 ; i<n ; i++) {
177 for (int j=0 ; j<n ; j++)
178 vect[j] = 0 ;
179 vect[i] = 1 ;
180 d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
181 vect[i] += masse*masse*alpha*alpha ;
182 multx_1d (n, &vect, R_CHEB) ;
183 multx_1d (n, &vect, R_CHEB) ;
184 for (int j=0 ; j<n ; j++)
185 dd.set(j, i) += vect[j] ;
186 }
187
188 for (int i=0 ; i<n ; i++) {
189 for (int j=0 ; j<n ; j++)
190 vect[j] = 0 ;
191 vect[i] = 1 ;
192 sxdsdx_1d (n, &vect, R_CHEB) ;
193 for (int j=0 ; j<n ; j++)
194 xd.set(j, i) = vect[j]*echelle ;
195 }
196
197 for (int i=0 ; i<n ; i++) {
198 for (int j=0 ; j<n ; j++)
199 vect[j] = 0 ;
200 vect[i] = 1 ;
201 sxdsdx_1d (n, &vect, R_CHEB) ;
202 multx_1d (n, &vect, R_CHEB) ;
203 for (int j=0 ; j<n ; j++)
204 xd.set(j, i) += vect[j] ;
205 }
206
207 for (int i=0 ; i<n ; i++) {
208 for (int j=0 ; j<n ; j++)
209 vect[j] = 0 ;
210 vect[i] = 1 ;
211 sx2_1d (n, &vect, R_CHEB) ;
212 for (int j=0 ; j<n ; j++)
213 xx.set(j, i) = vect[j] ;
214 }
215
216 delete [] vect ;
217
218 Matrice res(n, n) ;
219 res = dd+2*xd - lq*(lq+1)*xx;
220
221 return res ;
222}
223
224
225 //--------------------------
226 //- La routine a appeler ---
227 //----------------------------
228
229Matrice helmholtz_plus_mat(int n, int lq, double alpha, double beta, double masse,
230 int base_r)
231{
232
233 // Routines de derivation
234 static Matrice (*helmholtz_plus_mat[MAX_BASE])(int, int, double, double, double);
235 static int nap = 0 ;
236
237 // Premier appel
238 if (nap==0) {
239 nap = 1 ;
240 for (int i=0 ; i<MAX_BASE ; i++) {
241 helmholtz_plus_mat[i] = _helmholtz_plus_mat_pas_prevu ;
242 }
243 // Les routines existantes
244 helmholtz_plus_mat[R_CHEB >> TRA_R] = _helmholtz_plus_mat_r_cheb ;
245 helmholtz_plus_mat[R_CHEBP >> TRA_R] = _helmholtz_plus_mat_r_chebp ;
246 }
247
248 Matrice res(helmholtz_plus_mat[base_r](n, lq, alpha, beta, masse)) ;
249 return res ;
250}
251
252}
#define MAX_BASE
Nombre max. de bases differentes.
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
#define R_CHEB
base de Chebychev ordinaire (fin)
#define R_CHEBP
base de Cheb. paire (rare) seulement
Lorene prototypes.
Definition app_hor.h:64