LORENE
val_solp.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 val_solp_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/val_solp.C,v 1.6 2014/10/13 08:53:31 j_novak Exp $" ;
24
25/*
26 * $Id: val_solp.C,v 1.6 2014/10/13 08:53:31 j_novak Exp $
27 * $Log: val_solp.C,v $
28 * Revision 1.6 2014/10/13 08:53:31 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.5 2014/10/06 15:16:11 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.4 2008/02/18 13:53:45 j_novak
35 * Removal of special indentation instructions.
36 *
37 * Revision 1.3 2004/08/24 09:14:44 p_grandclement
38 * Addition of some new operators, like Poisson in 2d... It now requieres the
39 * GSL library to work.
40 *
41 * Also, the way a variable change is stored by a Param_elliptic is changed and
42 * no longer uses Change_var but rather 2 Scalars. The codes using that feature
43 * will requiere some modification. (It should concern only the ones about monopoles)
44 *
45 * Revision 1.2 2003/12/11 15:37:09 p_grandclement
46 * sqrt(2) ----> sqrt(double(2))
47 *
48 * Revision 1.1 2003/12/11 14:48:49 p_grandclement
49 * Addition of ALL (and that is a lot !) the files needed for the general elliptic solver ... UNDER DEVELOPEMENT...
50 *
51 *
52 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/val_solp.C,v 1.6 2014/10/13 08:53:31 j_novak Exp $
53 *
54 */
55
56//fichiers includes
57#include <cstdio>
58#include <cstdlib>
59#include <cmath>
60
61#include "proto.h"
62#include "matrice.h"
63#include "type_parite.h"
64
65
66 //------------------------------------
67 // Routine pour les cas non prevus --
68 //------------------------------------
69namespace Lorene {
70Tbl _val_solp_pas_prevu (const Tbl&, double) {
71
72 cout << " Base_r unknown in val_solp."<< endl ;
73 abort() ;
74 exit(-1) ;
75 Tbl res(1) ;
76 return res;
77}
78
79
80 //-------------------
81 //-- R_CHEB ------
82 //-------------------
83
84Tbl _val_solp_r_cheb (const Tbl& sp, double alpha) {
85
86 int nr = sp.get_dim(0) ;
87 Tbl res(4) ;
88 res.annule_hard() ;
89
90 // Solution en + 1
91 for (int i=0 ; i<nr ; i++)
92 res.set(0) += sp(i) ;
93
94 // Solution en -1 :
95 for (int i=0 ; i<nr ; i++)
96 if (i%2 == 0)
97 res.set(1) += sp(i) ;
98 else
99 res.set(1) -= sp(i) ;
100
101 // Derivee en +1 :
102 for (int i=0 ; i<nr ; i++)
103 res.set(2) += sp(i)*i*i/alpha ;
104
105 // Derivee en -1 :
106 for (int i=0 ; i<nr ; i++)
107 if (i%2 == 0)
108 res.set(3) -= sp(i)*i*i/alpha ;
109 else
110 res.set(3) += sp(i)*i*i/alpha ;
111
112 res /= sqrt(double(2)) ;
113 return res ;
114}
115
116 //-------------------
117 //-- R_CHEBP ------
118 //-------------------
119
120Tbl _val_solp_r_chebp (const Tbl& sp, double alpha) {
121
122 int nr = sp.get_dim(0) ;
123 Tbl res(4) ;
124 res.annule_hard() ;
125
126 // Solution en +1 :
127 for (int i=0 ; i<nr ; i++)
128 res.set(0) += sp(i) ;
129
130 // Solution en 0 (a priori pas trop utilise)
131 for (int i=0 ; i<nr ; i++)
132 if (i%2==0)
133 res.set(1) += sp(i) ;
134 else
135 res.set(1) -= sp(i) ;
136
137 // Derivee en +1 :
138 for (int i=0 ; i<nr ; i++)
139 res.set(2) += sp(i)*(2*i)*(2*i)/alpha ;
140
141 // Derivee en 0
142 res.set(3) = 0 ;
143
144 res /= sqrt(double(2)) ;
145 return res ;
146}
147
148
149 //-------------------
150 //-- R_CHEBI -----
151 //-------------------
152
153Tbl _val_solp_r_chebi (const Tbl& sp, double alpha) {
154
155 int nr = sp.get_dim(0) ;
156 Tbl res(4) ;
157 res.annule_hard() ;
158
159 // Solution en +1 :
160 for (int i=0 ; i<nr ; i++)
161 res.set(0) += sp(i) ;
162
163 // Solution en 0 :
164 res.set(1) = 0 ;
165
166 // Derivee en +1 :
167 for (int i=0 ; i<nr ; i++)
168 res.set(2) += sp(i)*(2*i+1)*(2*i+1)/alpha ;
169
170 // Derivee en 0 :
171 for (int i=0 ; i<nr ; i++)
172 if (i%2==0)
173 res.set(3) += (2*i+1)*sp(i) ;
174 else
175 res.set(3) -= (2*i+1)*sp(i) ;
176
177 res /= sqrt(double(2)) ;
178 return res ;
179}
180
181
182
183 //-------------------
184 //-- R_CHEBU -----
185 //-------------------
186
187Tbl _val_solp_r_chebu (const Tbl& sp, double alpha) {
188
189 int nr = sp.get_dim(0) ;
190 Tbl res(4) ;
191 res.annule_hard() ;
192
193 // Solution en + 1
194 for (int i=0 ; i<nr ; i++)
195 res.set(0) += sp(i) ;
196
197 // Solution en -1 :
198 for (int i=0 ; i<nr ; i++)
199 if (i%2==0)
200 res.set(1) += sp(i) ;
201 else
202 res.set(1) -= sp(i) ;
203
204 // Derivee en +1 c'est zero ca !
205
206 // Derivee en -1 :
207 for (int i=0 ; i<nr ; i++)
208 if (i%2==0)
209 res.set(3) += 4.*alpha*i*i*sp(i) ;
210 else
211 res.set(3) -= 4.*alpha*i*i*sp(i) ;
212
213 res /= sqrt(double(2)) ;
214 return res ;
215}
216
217
218
219
220 //-------------------
221 //-- Fonction ----
222 //-------------------
223
224
225Tbl val_solp (const Tbl& sp, double alpha, int base_r) {
226
227 // Routines de derivation
228 static Tbl (*val_solp[MAX_BASE])(const Tbl&, double) ;
229 static int nap = 0 ;
230
231 // Premier appel
232 if (nap==0) {
233 nap = 1 ;
234 for (int i=0 ; i<MAX_BASE ; i++) {
235 val_solp[i] = _val_solp_pas_prevu ;
236 }
237 // Les routines existantes
238 val_solp[R_CHEB >> TRA_R] = _val_solp_r_cheb ;
239 val_solp[R_CHEBU >> TRA_R] = _val_solp_r_chebu ;
240 val_solp[R_CHEBP >> TRA_R] = _val_solp_r_chebp ;
241 val_solp[R_CHEBI >> TRA_R] = _val_solp_r_chebi ;
242 }
243
244 Tbl res(val_solp[base_r](sp, alpha)) ;
245 return res ;
246}
247}
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition tbl.h:403
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
#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