LORENE
solh_helmholtz_plus.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 solh_helmholtz_plusC[] = "$Header $" ;
24
25
26/*
27 * $Id: solh_helmholtz_plus.C,v 1.6 2014/10/13 08:53:31 j_novak Exp $
28 * $Log: solh_helmholtz_plus.C,v $
29 * Revision 1.6 2014/10/13 08:53:31 j_novak
30 * Lorene classes and functions now belong to the namespace Lorene.
31 *
32 * Revision 1.5 2014/10/06 15:16:10 j_novak
33 * Modified #include directives to use c++ syntax.
34 *
35 * Revision 1.4 2008/02/18 13:53:43 j_novak
36 * Removal of special indentation instructions.
37 *
38 * Revision 1.3 2007/05/06 10:48:12 p_grandclement
39 * Modification of a few operators for the vorton project
40 *
41 * Revision 1.2 2004/01/15 09:15:37 p_grandclement
42 * Modification and addition of the Helmholtz operators
43 *
44 * Revision 1.1 2003/12/11 14:48:49 p_grandclement
45 * Addition of ALL (and that is a lot !) the files needed for the general elliptic solver ... UNDER DEVELOPEMENT...
46 *
47 *
48 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/solh_helmholtz_plus.C,v 1.6 2014/10/13 08:53:31 j_novak Exp $
49 *
50 */
51
52//fichiers includes
53#include <cstdio>
54#include <cstdlib>
55#include <cmath>
56#include <gsl/gsl_sf_bessel.h>
57
58#include "proto.h"
59#include "matrice.h"
60#include "type_parite.h"
61
62 //------------------------------------
63 // Routine pour les cas non prevus --
64 //------------------------------------
65namespace Lorene {
66Tbl _solh_helmholtz_plus_pas_prevu (int, int, double, double, double) {
67
68 cout << "Homogeneous solution not implemented in hemlholtz_plus : "<< endl ;
69 abort() ;
70 exit(-1) ;
71 Tbl res(1) ;
72 return res;
73}
74
75
76
77 //-------------------
78 //-- R_CHEB ------
79 //-------------------
80
81Tbl _solh_helmholtz_plus_r_cheb (int n, int lq, double alpha, double beta,
82 double masse) {
83
84 assert (masse > 0) ;
85
86 Tbl res(2,n) ;
87 res.set_etat_qcq() ;
88 double* coloc = new double[n] ;
89
90 int * deg = new int[3] ;
91 deg[0] = 1 ;
92 deg[1] = 1 ;
93 deg[2] = n ;
94
95 // First SH
96 for (int i=0 ; i<n ; i++){
97 double air = alpha*(-cos(M_PI*i/(n-1))) + beta ;
98 coloc[i] = -gsl_sf_bessel_jl(lq, masse*air) ;
99 }
100
101 cfrcheb(deg, deg, coloc, deg, coloc) ;
102 for (int i=0 ; i<n ;i++)
103 res.set(0,i) = coloc[i] ;
104
105 // Second SH
106 for (int i=0 ; i<n ; i++){
107 double air = alpha*(-cos(M_PI*i/(n-1))) + beta ;
108 coloc[i] = -gsl_sf_bessel_yl (lq, masse*air) ;
109 }
110
111 cfrcheb(deg, deg, coloc, deg, coloc) ;
112 for (int i=0 ; i<n ;i++)
113 res.set(1,i) = coloc[i] ;
114
115 delete [] deg ;
116 delete [] coloc ;
117 return res ;
118}
119
120
121 //-------------------
122 //-- R_CHEBP -----
123 //-------------------
124
125Tbl _solh_helmholtz_plus_r_chebp (int n, int lq, double alpha, double,
126 double masse) {
127
128 assert (masse > 0) ;
129
130 Tbl res(n) ;
131 res.set_etat_qcq() ;
132 double* coloc = new double[n] ;
133
134 int * deg = new int[3] ;
135 deg[0] = 1 ;
136 deg[1] = 1 ;
137 deg[2] = n ;
138
139 // SH
140 for (int i=0 ; i<n ; i++){
141 double air = alpha*sin(M_PI*i/2/(n-1)) ;
142 coloc[i] = -gsl_sf_bessel_jl(lq, masse*air) ;
143 }
144
145 cfrchebp(deg, deg, coloc, deg, coloc) ;
146 for (int i=0 ; i<n ;i++)
147 res.set(i) = coloc[i] ;
148
149 delete [] deg ;
150 delete [] coloc ;
151 return res ;
152}
153
154
155 //-------------------
156 //-- Fonction ----
157 //-------------------
158
159
160Tbl solh_helmholtz_plus (int n, int lq, double alpha, double beta,
161 double masse, int base_r) {
162
163 // Routines de derivation
164 static Tbl (*solh_helmholtz_plus[MAX_BASE])(int, int, double, double, double) ;
165 static int nap = 0 ;
166
167 // Premier appel
168 if (nap==0) {
169 nap = 1 ;
170 for (int i=0 ; i<MAX_BASE ; i++) {
171 solh_helmholtz_plus[i] = _solh_helmholtz_plus_pas_prevu ;
172 }
173 // Les routines existantes
174 solh_helmholtz_plus[R_CHEB >> TRA_R] = _solh_helmholtz_plus_r_cheb ;
175 solh_helmholtz_plus[R_CHEBP >> TRA_R] = _solh_helmholtz_plus_r_chebp ;
176 }
177
178 Tbl res(solh_helmholtz_plus[base_r](n, lq, alpha, beta, masse)) ;
179 return res ;
180}
181}
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:69
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
#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