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
23char 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 //------------------------------------
61namespace Lorene {
62Tbl _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
76Tbl _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
103Tbl _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
123Tbl _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
144Tbl _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
166Tbl 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}
Tbl & set(int l)
Read/write of the value in a given domain.
Definition cmp.h:724
Basic array class.
Definition tbl.h:161
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
#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
Tbl _val_solh_r_chebu(int l, double alpha, double)
Definition val_solh.C:144