LORENE
comb_lin_helmholtz_minus.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 comb_lin_helmholtz_minusC[] = "$Header $" ;
24
25/*
26 * $Id: comb_lin_helmholtz_minus.C,v 1.7 2014/10/13 08:53:28 j_novak Exp $
27 * $Log: comb_lin_helmholtz_minus.C,v $
28 * Revision 1.7 2014/10/13 08:53:28 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.6 2014/10/06 15:16:08 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.5 2008/07/08 11:45:28 p_grandclement
35 * Add helmholtz_minus in the nucleus
36 *
37 * Revision 1.4 2008/02/18 13:53:42 j_novak
38 * Removal of special indentation instructions.
39 *
40 * Revision 1.3 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.2 2004/01/15 09:15:37 p_grandclement
49 * Modification and addition of the Helmholtz operators
50 *
51 * Revision 1.1 2003/12/11 14:48:49 p_grandclement
52 * Addition of ALL (and that is a lot !) the files needed for the general elliptic solver ... UNDER DEVELOPEMENT...
53 *
54 *
55 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/comb_lin_helmholtz_minus.C,v 1.7 2014/10/13 08:53:28 j_novak Exp $
56 *
57 */
58
59//fichiers includes
60#include <cstdio>
61#include <cstdlib>
62#include <cmath>
63
64#include "matrice.h"
65#include "type_parite.h"
66#include "proto.h"
67
68
69// Version Matrice --> Matrice
70namespace Lorene {
71Matrice _cl_helmholtz_minus_pas_prevu (const Matrice& so) {
72 cout << "CL Helmholtz minus not implemented" << endl ;
73 abort() ;
74 exit(-1) ;
75 return so;
76}
77
78
79
80 //-------------------
81 //-- R_CHEB ------
82 //-------------------
83
84Matrice _cl_helmholtz_minus_r_cheb (const Matrice& source) {
85
86 int n = source.get_dim(0) ;
87 assert (n==source.get_dim(1)) ;
88
89 Matrice barre(source) ;
90 int dirac = 1 ;
91 for (int i=0 ; i<n-2 ; i++) {
92 for (int j=0 ; j<n ; j++)
93 barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j))
94 /(i+1) ;
95 if (i==0) dirac = 0 ;
96 }
97
98 Matrice res(barre) ;
99 for (int i=0 ; i<n-4 ; i++)
100 for (int j=0 ; j<n ; j++)
101 res.set(i, j) = barre(i, j)-barre(i+2, j) ;
102
103 return res ;
104}
105
106
107 //-------------------
108 //-- R_CHEBU ------
109 //-------------------
110
111Matrice _cl_helmholtz_minus_r_chebu (const Matrice& source) {
112
113 int n = source.get_dim(0) ;
114 assert (n==source.get_dim(1)) ;
115
116 Matrice barre(source) ;
117 int dirac = 1 ;
118 for (int i=0 ; i<n-2 ; i++) {
119 for (int j=0 ; j<n ; j++)
120 barre.set(i, j) = ((1+dirac)*source(i, j)-source(i+2, j)) ;
121 if (i==0) dirac = 0 ;
122 }
123
124 Matrice tilde(barre) ;
125 for (int i=0 ; i<n-4 ; i++)
126 for (int j=0 ; j<n ; j++)
127 tilde.set(i, j) = (barre(i, j)-barre(i+2, j)) ;
128
129 Matrice hat(tilde) ;
130 for (int i=0 ; i<n-4 ; i++)
131 for (int j=0 ; j<n ; j++)
132 hat.set(i, j) = (tilde(i, j)+tilde(i+1, j)) ;
133
134 Matrice res(hat) ;
135 for (int i=0 ; i<n-4 ; i++)
136 for (int j=0 ; j<n ; j++)
137 res.set(i, j) = hat(i, j)-hat(i+1, j) ;
138
139 return res ;
140}
141
142 //-------------------
143 //-- R_CHEBP -----
144 //-------------------
145
146
147Matrice _cl_helmholtz_minus_r_chebp (const Matrice &source) {
148 int n = source.get_dim(0) ;
149 assert (n==source.get_dim(1)) ;
150
151 Matrice barre(source) ;
152
153 int dirac = 1 ;
154 for (int i=0 ; i<n-2 ; i++) {
155 for (int j=0 ; j<n ; j++)
156 barre.set(i, j) = (1+dirac)*source(i, j)-source(i+2, j) ;
157 if (i==0) dirac = 0 ;
158 }
159
160 Matrice tilde(barre) ;
161 for (int i=0 ; i<n-4 ; i++)
162 for (int j=0 ; j<n ; j++)
163 tilde.set(i, j) = barre(i, j)-barre(i+2, j) ;
164
165 Matrice res(tilde) ;
166 for (int i=0 ; i<n-4 ; i++)
167 for (int j=0 ; j<n ; j++)
168 res.set(i, j) = tilde(i, j)-tilde(i+1, j) ;
169 return res ;
170}
171
172 //-------------------
173 //-- R_CHEBI -----
174 //-------------------
175
176
177Matrice _cl_helmholtz_minus_r_chebi (const Matrice &source) {
178 int n = source.get_dim(0) ;
179 assert (n==source.get_dim(1)) ;
180
181 Matrice barre(source) ;
182
183 for (int i=0 ; i<n-2 ; i++)
184 for (int j=0 ; j<n ; j++)
185 barre.set(i, j) = source(i, j)-source(i+2, j) ;
186
187 Matrice tilde(barre) ;
188 for (int i=0 ; i<n-4 ; i++)
189 for (int j=0 ; j<n ; j++)
190 tilde.set(i, j) = barre(i, j)-barre(i+2, j) ;
191
192 Matrice res(tilde) ;
193 for (int i=0 ; i<n-4 ; i++)
194 for (int j=0 ; j<n ; j++)
195 res.set(i, j) = tilde(i, j)-tilde(i+1, j) ;
196 return res ;
197}
198
199 //-------------------------
200 //- La routine a appeler ---
201 //---------------------------
202
203Matrice cl_helmholtz_minus (const Matrice &source, int base_r) {
204
205 // Routines de derivation
206 static Matrice (*cl_helmholtz_minus[MAX_BASE]) (const Matrice &) ;
207 static int nap = 0 ;
208
209 // Premier appel
210 if (nap==0) {
211 nap = 1 ;
212 for (int i=0 ; i<MAX_BASE ; i++) {
213 cl_helmholtz_minus[i] = _cl_helmholtz_minus_pas_prevu ;
214 }
215 // Les routines existantes
216 cl_helmholtz_minus[R_CHEB >> TRA_R] = _cl_helmholtz_minus_r_cheb ;
217 cl_helmholtz_minus[R_CHEBU >> TRA_R] = _cl_helmholtz_minus_r_chebu ;
218 cl_helmholtz_minus[R_CHEBP >> TRA_R] = _cl_helmholtz_minus_r_chebp ;
219 cl_helmholtz_minus[R_CHEBI >> TRA_R] = _cl_helmholtz_minus_r_chebi ;
220 }
221
222 Matrice res(cl_helmholtz_minus[base_r](source)) ;
223 return res ;
224}
225
226
227//************************ TBL Versions *************************************
228
229
230
231
232Tbl _cl_helmholtz_minus_pas_prevu (const Tbl &so) {
233
234 cout << "Linear combination for Helmholtz minus not implemented..." << endl ;
235 abort() ;
236 exit(-1) ;
237 return so;
238}
239
240 //-------------------
241 //-- R_CHEB -------
242 //--------------------
243Tbl _cl_helmholtz_minus_r_cheb (const Tbl& source) {
244
245 int n = source.get_dim(0) ;
246
247 Tbl barre(source) ;
248 int dirac = 1 ;
249 for (int i=0 ; i<n-2 ; i++) {
250 barre.set(i) = ((1+dirac)*source(i)-source(i+2))
251 /(i+1) ;
252 if (i==0) dirac = 0 ;
253 }
254
255 Tbl res(barre) ;
256 for (int i=0 ; i<n-4 ; i++)
257 res.set(i) = barre(i)-barre(i+2) ;
258
259 return res ;
260}
261
262
263 //------------------
264 //-- R_CHEBU -------
265 //--------------------
266
267Tbl _cl_helmholtz_minus_r_chebu (const Tbl& source) {
268
269 int n = source.get_dim(0) ;
270
271 Tbl barre(source) ;
272 int dirac = 1 ;
273 for (int i=0 ; i<n-2 ; i++) {
274 barre.set(i) = ((1+dirac)*source(i)-source(i+2)) ;
275 if (i==0) dirac = 0 ;
276 }
277
278 Tbl tilde(barre) ;
279 for (int i=0 ; i<n-4 ; i++)
280 tilde.set(i) = (barre(i)-barre(i+2)) ;
281
282 Tbl hat(tilde) ;
283 for (int i=0 ; i<n-4 ; i++)
284 hat.set(i) = (tilde(i)+tilde(i+1)) ;
285
286 Tbl res(hat) ;
287 for (int i=0 ; i<n-4 ; i++)
288 res.set(i) = hat(i)-hat(i+1) ;
289
290 return res ;
291}
292
293 //-------------------
294 //-- R_CHEBP -----
295 //-------------------
296
297
298Tbl _cl_helmholtz_minus_r_chebp (const Tbl &source) {
299 int n = source.get_dim(0) ;
300 Tbl barre(source) ;
301
302 int dirac = 1 ;
303 for (int i=0 ; i<n-2 ; i++) {
304 barre.set(i) = (1+dirac)*source(i)-source(i+2) ;
305 if (i==0) dirac = 0 ;
306 }
307
308 Tbl tilde(barre) ;
309 for (int i=0 ; i<n-4 ; i++)
310 tilde.set(i) = barre(i)-barre(i+2) ;
311
312 Tbl res(tilde) ;
313 for (int i=0 ; i<n-4 ; i++)
314 res.set(i) = tilde(i)-tilde(i+1) ;
315 return res ;
316}
317
318 //-------------------
319 //-- R_CHEBI -----
320 //-------------------
321
322
323Tbl _cl_helmholtz_minus_r_chebi (const Tbl &source) {
324 int n = source.get_dim(0) ;
325 Tbl barre(source) ;
326
327 for (int i=0 ; i<n-2 ; i++)
328 barre.set(i) = source(i)-source(i+2) ;
329
330 Tbl tilde(barre) ;
331 for (int i=0 ; i<n-4 ; i++)
332 tilde.set(i) = barre(i)-barre(i+2) ;
333
334 Tbl res(tilde) ;
335 for (int i=0 ; i<n-4 ; i++)
336 res.set(i) = tilde(i)-tilde(i+1) ;
337 return res ;
338}
339 //----------------------------
340 //- Routine a appeler ---
341 //------------------------------
342
343Tbl cl_helmholtz_minus (const Tbl &source, int base_r) {
344
345 // Routines de derivation
346 static Tbl (*cl_helmholtz_minus[MAX_BASE])(const Tbl &) ;
347 static int nap = 0 ;
348
349 // Premier appel
350 if (nap==0) {
351 nap = 1 ;
352 for (int i=0 ; i<MAX_BASE ; i++) {
353 cl_helmholtz_minus[i] = _cl_helmholtz_minus_pas_prevu ;
354 }
355 // Les routines existantes
356 cl_helmholtz_minus[R_CHEB >> TRA_R] = _cl_helmholtz_minus_r_cheb ;
357 cl_helmholtz_minus[R_CHEBU >> TRA_R] = _cl_helmholtz_minus_r_chebu ;
358 cl_helmholtz_minus[R_CHEBP >> TRA_R] = _cl_helmholtz_minus_r_chebp ;
359 cl_helmholtz_minus[R_CHEBI >> TRA_R] = _cl_helmholtz_minus_r_chebi ;
360 }
361
362 Tbl res(cl_helmholtz_minus[base_r](source)) ;
363 return res ;
364}
365}
int get_dim(int i) const
Returns the dimension of the matrix.
Definition matrice.C:260
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#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