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