LORENE
val_solh.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_solh_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/val_solh.C,v 1.5 2014/10/13 08:53:31 j_novak Exp $" ;
24 
25 /*
26  * $Id: val_solh.C,v 1.5 2014/10/13 08:53:31 j_novak Exp $
27  * $Log: val_solh.C,v $
28  * Revision 1.5 2014/10/13 08:53:31 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.4 2014/10/06 15:16:11 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.3 2008/02/18 13:53:45 j_novak
35  * Removal of special indentation instructions.
36  *
37  * Revision 1.2 2003/12/11 15:37:09 p_grandclement
38  * sqrt(2) ----> sqrt(double(2))
39  *
40  * Revision 1.1 2003/12/11 14:48:49 p_grandclement
41  * Addition of ALL (and that is a lot !) the files needed for the general elliptic solver ... UNDER DEVELOPEMENT...
42  *
43  *
44  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/val_solh.C,v 1.5 2014/10/13 08:53:31 j_novak Exp $
45  *
46  */
47 
48 //fichiers includes
49 #include <cstdio>
50 #include <cstdlib>
51 #include <cmath>
52 
53 #include "proto.h"
54 #include "matrice.h"
55 #include "type_parite.h"
56 
57 
58  //------------------------------------
59  // Routine pour les cas non prevus --
60  //------------------------------------
61 namespace Lorene {
62 Tbl _val_solh_pas_prevu (int, double, double) {
63 
64  cout << " Solution homogene pas prevue ..... : "<< endl ;
65  abort() ;
66  exit(-1) ;
67  Tbl res(1) ;
68  return res;
69 }
70 
71 
72  //-------------------
73  //-- R_CHEB ------
74  //-------------------
75 
76 Tbl _val_solh_r_cheb (int l, double alpha, double beta) {
77 
78  double echelle = beta/alpha ;
79 
80  Tbl res(2, 4) ;
81  res.set_etat_qcq() ;
82 
83  // Solution 1 : (x+echelle)^l
84  res.set(0,0) = pow(1.+echelle, l) ;
85  res.set(0,1) = pow(-1.+echelle, l) ;
86  res.set(0,2) = pow(1.+echelle, l-1)*l/alpha ;
87  res.set(0,3) = pow(-1.+echelle, l-1)*l/alpha ;
88 
89  // Solution 2 : 1./(x+echelle)^(l+1)
90  res.set(1,0) = 1./pow(1.+echelle, l+1) ;
91  res.set(1,1) = 1./pow(-1.+echelle, l+1) ;
92  res.set(1,2) = -1./pow(1.+echelle, l+2)*(l+1)/alpha ;
93  res.set(1,3) = -1./pow(-1.+echelle, l+2)*(l+1)/alpha ;
94 
95  res /= sqrt(double(2)) ;
96  return res ;
97 }
98 
99  //-------------------
100  //-- R_CHEBP ------
101  //-------------------
102 
103 Tbl _val_solh_r_chebp (int l, double alpha, double) {
104 
105  Tbl res(4) ;
106  res.set_etat_qcq() ;
107 
108  // Solution : x^l
109  res.set(0) = 1. ;
110  res.set(1) = (l==0) ? 1. : 0. ;
111  res.set(2) = 1./alpha*l ;
112  res.set(3) = (l==1) ? 1 : 0 ;
113 
114  res /= sqrt(double(2)) ;
115  return res ;
116 }
117 
118 
119  //-------------------
120  //-- R_CHEBI -----
121  //-------------------
122 
123 Tbl _val_solh_r_chebi (int l, double alpha, double) {
124 
125  Tbl res(4) ;
126  res.set_etat_qcq() ;
127 
128  // Solution : x^l
129  res.set(0) = 1. ;
130  res.set(1) = (l==0) ? 1. : 0 ;
131  res.set(2) = 1./alpha*l ;
132  res.set(3) = (l==1) ? 1 : 0. ;
133 
134  res /= sqrt(double(2)) ;
135  return res ;
136 }
137 
138 
139 
140  //-------------------
141  //-- R_CHEBU -----
142  //-------------------
143 
144 Tbl _val_solh_r_chebu (int l, double alpha, double) {
145 
146  Tbl res(4) ;
147  res.set_etat_qcq() ;
148 
149  // Solution : 1/(x-1)^(l+1)
150  res.set(0) = 0 ;
151  res.set(1) = pow(-2., l+1)/sqrt(double(2)) ;
152  res.set(2) = 0. ;
153  res.set(3) = -alpha*(l+1)*pow(-2., l+2)/sqrt(double(2)) ;
154 
155  return res ;
156 }
157 
158 
159 
160 
161  //-------------------
162  //-- Fonction ----
163  //-------------------
164 
165 
166 Tbl val_solh(int l, double alpha, double beta, int base_r) {
167 
168  // Routines de derivation
169  static Tbl (*val_solh[MAX_BASE])(int, double, double) ;
170  static int nap = 0 ;
171 
172  // Premier appel
173  if (nap==0) {
174  nap = 1 ;
175  for (int i=0 ; i<MAX_BASE ; i++) {
176  val_solh[i] = _val_solh_pas_prevu ;
177  }
178  // Les routines existantes
179  val_solh[R_CHEB >> TRA_R] = _val_solh_r_cheb ;
180  val_solh[R_CHEBU >> TRA_R] = _val_solh_r_chebu ;
181  val_solh[R_CHEBP >> TRA_R] = _val_solh_r_chebp ;
182  val_solh[R_CHEBI >> TRA_R] = _val_solh_r_chebi ;
183  }
184 
185  Tbl res(val_solh[base_r](l, alpha, beta)) ;
186  return res ;
187 }
188 }
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::_val_solh_r_chebu
Tbl _val_solh_r_chebu(int l, double alpha, double)
Definition: val_solh.C:144
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
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::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::Cmp::set
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:724
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