LORENE
ope_sec_order_solh.C
1/*
2 * Copyright (c) 2003 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 version 2
8 * as published by the Free Software Foundation.
9 *
10 * LORENE is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with LORENE; if not, write to the Free Software
17 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 *
19 */
20
21char ope_sec_order_solh_C[] = "$Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_sec_order/ope_sec_order_solh.C,v 1.4 2014/10/13 08:53:36 j_novak Exp $" ;
22
23/*
24 * $Id: ope_sec_order_solh.C,v 1.4 2014/10/13 08:53:36 j_novak Exp $
25 * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_sec_order/ope_sec_order_solh.C,v 1.4 2014/10/13 08:53:36 j_novak Exp $
26 *
27 */
28#include <cmath>
29#include <cstdlib>
30
31#include "proto.h"
32#include "ope_elementary.h"
33
34 //------------------------------------
35 // Routine pour les cas non prevus --
36 //------------------------------------
37namespace Lorene {
38Tbl _solh_sec_order_pas_prevu (int, double, double,double,double,double,Tbl&) {
39
40 cout << "Homogeneous solution not implemented in Sec_order : "<< endl ;
41 abort() ;
42 exit(-1) ;
43 Tbl res(1) ;
44 return res;
45}
46
47
48 //-------------------
49 //-- R_CHEB ------
50 //-------------------
51
52Tbl _solh_sec_order_r_cheb (int n, double alpha, double beta,
53 double a, double b, double c, Tbl& val_lim) {
54
55 // Stuff to compute the coefs...
56 Tbl res(2,n) ;
57 res.set_etat_qcq() ;
58 double* coloc = new double[n] ;
59
60 int * deg = new int[3] ;
61 deg[0] = 1 ;
62 deg[1] = 1 ;
63 deg[2] = n ;
64
65 // Array on the radius
66 double* sigma = new double[n] ;
67 for (int i=0 ; i<n ; i++)
68 sigma[i] = alpha*(-cos(M_PI*i/(n-1))) + beta ;
69
70 double sigma_minus = beta-alpha ;
71 double sigma_plus = beta+alpha ;
72
73 // Determinant :
74 double delta = b*b-4*a*c ;
75 int signe = 0 ;
76 if (delta > 0)
77 signe = + 1 ;
78 if (delta < 0)
79 signe = -1 ;
80 if (fabs(delta) < 1e-14)
81 signe = 0 ;
82
83 switch (signe) {
84
85 case 1: {
86 // Two real solutions
87 double lambda_one = (-b + sqrt(delta)) / 2./a ;
88 double lambda_two = (-b - sqrt(delta)) / 2./a ;
89
90 // First SH
91 for (int i=0 ; i<n ; i++)
92 coloc[i] = exp(sigma[i]*lambda_one) ;
93 cfrcheb(deg, deg, coloc, deg, coloc) ;
94 for (int i=0 ; i<n ;i++)
95 res.set(0,i) = coloc[i] ;
96
97 // Second SH
98 for (int i=0 ; i<n ; i++)
99 coloc[i] = exp(sigma[i]*lambda_two) ;
100 cfrcheb(deg, deg, coloc, deg, coloc) ;
101 for (int i=0 ; i<n ;i++)
102 res.set(1,i) = coloc[i] ;
103
104 // Limit on the boundaries :
105 val_lim.set(0,0) = exp(sigma_minus*lambda_one) ;
106 val_lim.set(0,1) = lambda_one*exp(sigma_minus*lambda_one)
107 /exp(sigma_minus) ;
108 val_lim.set(0,2) = exp(sigma_plus*lambda_one) ;
109 val_lim.set(0,3) = lambda_one*exp(sigma_plus*lambda_one)
110 /exp(sigma_plus) ;
111
112 val_lim.set(1,0) = exp(sigma_minus*lambda_two) ;
113 val_lim.set(1,1) = lambda_two*exp(sigma_minus*lambda_two)/
114 exp(sigma_minus);
115 val_lim.set(1,2) = exp(sigma_plus*lambda_two) ;
116 val_lim.set(1,3) = lambda_two*exp(sigma_plus*lambda_two)
117 /exp(sigma_plus) ;
118 val_lim /= sqrt(double(2)) ;
119 break ;
120 }
121 case 0: {
122 // Only one solution :
123 double lambda = -b/2./a ;
124 // First SH
125 for (int i=0 ; i<n ; i++)
126 coloc[i] = exp(sigma[i]*lambda) ;
127 cfrcheb(deg, deg, coloc, deg, coloc) ;
128 for (int i=0 ; i<n ;i++)
129 res.set(0,i) = coloc[i] ;
130
131 // Second SH
132 for (int i=0 ; i<n ; i++)
133 coloc[i] = sigma[i]*exp(sigma[i]*lambda) ;
134 cfrcheb(deg, deg, coloc, deg, coloc) ;
135 for (int i=0 ; i<n ;i++)
136 res.set(1,i) = coloc[i] ;
137
138 // Limit on the boundaries :
139 val_lim.set(0,0) = exp(sigma_minus*lambda) ;
140 val_lim.set(0,1) = lambda*exp(sigma_minus*lambda)/exp(sigma_minus) ;
141 val_lim.set(0,2) = exp(sigma_plus*lambda) ;
142 val_lim.set(0,3) = lambda*exp(sigma_plus*lambda)/exp(sigma_plus) ;
143
144 val_lim.set(1,0) = sigma_minus*exp(sigma_minus*lambda) ;
145 val_lim.set(1,1) = exp(sigma_minus*lambda)* (lambda*sigma_minus+1)
146 /exp(sigma_minus);
147 val_lim.set(1,2) = sigma_plus*exp(sigma_plus*lambda) ;
148 val_lim.set(1,3) = exp(sigma_plus*lambda)* (lambda*sigma_plus+1)/
149 exp(sigma_plus) ;
150 val_lim /= sqrt(double(2)) ;
151 break ;
152 }
153 case -1:{
154 // Two imaginary solutions :
155 double real_part = -b/2./a ;
156 double imag_part = sqrt(-delta)/2./a ;
157
158 // First SH
159 for (int i=0 ; i<n ; i++)
160 coloc[i] = exp(sigma[i]*real_part)*cos(imag_part*sigma[i]) ;
161 cfrcheb(deg, deg, coloc, deg, coloc) ;
162 for (int i=0 ; i<n ;i++)
163 res.set(0,i) = coloc[i] ;
164
165 // Second SH
166 for (int i=0 ; i<n ; i++)
167 coloc[i] = exp(sigma[i]*real_part)*sin(imag_part*sigma[i]) ;
168 cfrcheb(deg, deg, coloc, deg, coloc) ;
169 for (int i=0 ; i<n ;i++)
170 res.set(1,i) = coloc[i] ;
171
172 // Limit on the boundaries :
173 val_lim.set(0,0) = exp(sigma_minus*real_part)*cos(imag_part*sigma_minus) ;
174 val_lim.set(0,1) = (real_part*cos(imag_part*sigma_minus) -
175 imag_part*sin(imag_part*sigma_minus)) * exp(real_part*sigma_minus)
176 /exp(sigma_minus);
177 val_lim.set(0,2) = exp(sigma_plus*real_part)*cos(imag_part*sigma_plus) ;
178 val_lim.set(0,3) = (real_part*cos(imag_part*sigma_plus) -
179 imag_part*sin(imag_part*sigma_plus)) * exp(real_part*sigma_plus)
180 /exp(sigma_plus) ;
181
182 val_lim.set(1,0) = exp(sigma_minus*real_part)*sin(imag_part*sigma_minus) ;
183 val_lim.set(1,1) = (real_part*sin(imag_part*sigma_minus) +
184 imag_part*cos(imag_part*sigma_minus)) * exp(real_part*sigma_minus)
185 /exp(sigma_minus);
186 val_lim.set(1,2) = exp(sigma_plus*real_part)*sin(imag_part*sigma_plus);
187 val_lim.set(1,3) = (real_part*sin(imag_part*sigma_plus) +
188 imag_part*cos(imag_part*sigma_plus)) * exp(real_part*sigma_plus)
189 /exp(sigma_plus) ;
190 val_lim /= sqrt(double(2)) ;
191 break ;
192 }
193 default:
194 cout << "What are you doing here ? Get out or I call the police !" << endl;
195 abort() ;
196 break ;
197 }
198
199 delete [] deg ;
200 delete [] coloc ;
201 delete [] sigma ;
202
203 return res ;
204}
205
206
208
209 // Routines de derivation
211 double, double, double, Tbl&) ;
212 static int nap = 0 ;
213
214 // Premier appel
215 if (nap==0) {
216 nap = 1 ;
217 for (int i=0 ; i<MAX_BASE ; i++) {
218 solh_sec_order[i] = _solh_sec_order_pas_prevu ;
219 }
220 // Les routines existantes
221 solh_sec_order[R_CHEB >> TRA_R] = _solh_sec_order_r_cheb ;
222 }
223
224 Tbl val_lim (2,4) ;
225 val_lim.set_etat_qcq() ;
227 val_lim)) ;
228
229 s_one_minus = val_lim(0,0) ;
230 ds_one_minus = val_lim(0,1) ;
231 s_one_plus = val_lim(0,2) ;
232 ds_one_plus = val_lim(0,3) ;
233
234 s_two_minus = val_lim(1,0) ;
235 ds_two_minus = val_lim(1,1) ;
236 s_two_plus = val_lim(1,2) ;
237 ds_two_plus = val_lim(1,3) ;
238
239 return res ;
240}
241}
Tbl & set(int l)
Read/write of the value in a given domain.
Definition cmp.h:724
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
double ds_two_minus
Value of the derivative of the second homogeneous solution at the inner boundary.
double s_two_plus
Value of the second homogeneous solution at the outer boundary.
double s_one_minus
Value of the first homogeneous solution at the inner boundary.
double beta
Parameter of the associated mapping.
double ds_one_plus
Value of the derivative of the first homogeneous solution at the outer boundary.
double ds_one_minus
Value of the derivative of the first homogeneous solution at the inner boundary.
double alpha
Parameter of the associated mapping.
double s_two_minus
Value of the second homogeneous solution at the inner boundary.
int base_r
Radial basis of decomposition.
double s_one_plus
Value of the first homogeneous solution at the outer boundary.
double ds_two_plus
Value of the derivative of the second homogeneous solution at the outer boundary.
int nr
Number of radial points.
double b_param
The parameter .
virtual Tbl get_solh() const
Computes the homogeneous solutions(s).
double a_param
The parameter .
double c_param
The parameter .
Basic array class.
Definition tbl.h:161
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:69
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
#define MAX_BASE
Nombre max. de bases differentes.
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
#define R_CHEB
base de Chebychev ordinaire (fin)
Lorene prototypes.
Definition app_hor.h:64