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 
23 char 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  //------------------------------------
69 namespace Lorene {
70 Tbl _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 
84 Tbl _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 
120 Tbl _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 
153 Tbl _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 
187 Tbl _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 
225 Tbl 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 }
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
Lorene
Lorene prototypes.
Definition: app_hor.h:64
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
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
Lorene::Tbl::get_dim
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:403
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
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