LORENE
ope_sec_order_r2_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_r2_solh_C[] = "$Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_sec_order_r2/ope_sec_order_r2_solh.C,v 1.4 2014/10/13 08:53:36 j_novak Exp $" ;
22
23/*
24 * $Id: ope_sec_order_r2_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_r2/ope_sec_order_r2_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_r2_pas_prevu (int, double, double,double,double,double,Tbl&) {
39
40 cout << "Homogeneous solution not implemented in Sec_order_r2 : "<< 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_r2_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* air = new double[n] ;
67 for (int i=0 ; i<n ; i++)
68 air[i] = alpha*(-cos(M_PI*i/(n-1))) + beta ;
69
70 double r_minus = beta-alpha ;
71 double r_plus = beta+alpha ;
72
73 // Determinant :
74 double delta = (b-a)*(b-a)-4*a*c ;
75 int signe ;
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 = ((a-b) + sqrt(delta)) / 2./a ;
88 double lambda_two = ((a-b) - sqrt(delta)) / 2./a ;
89
90 // First SH
91 for (int i=0 ; i<n ; i++)
92 coloc[i] = pow(air[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] = pow(air[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) = pow(r_minus, lambda_one) ;
106 val_lim.set(0,1) = lambda_one*pow(r_minus, lambda_one-1) ;
107 val_lim.set(0,2) = pow(r_plus, lambda_one) ;
108 val_lim.set(0,3) = lambda_one*pow(r_plus, lambda_one-1) ;
109
110 val_lim.set(1,0) = pow(r_minus, lambda_two) ;
111 val_lim.set(1,1) = lambda_two*pow(r_minus, lambda_two-1) ;
112 val_lim.set(1,2) = pow(r_plus, lambda_two) ;
113 val_lim.set(1,3) = lambda_two*pow(r_plus, lambda_two-1) ;
114 val_lim /= sqrt(double(2)) ;
115 break ;
116 }
117 case 0: {
118 // Only one solution :
119 double lambda = (a-b)/2./a ;
120 // First SH
121 for (int i=0 ; i<n ; i++)
122 coloc[i] = pow(air[i], lambda) ;
123 cfrcheb(deg, deg, coloc, deg, coloc) ;
124 for (int i=0 ; i<n ;i++)
125 res.set(0,i) = coloc[i] ;
126
127 // Second SH
128 for (int i=0 ; i<n ; i++)
129 coloc[i] = log(air[i])*pow(air[i], lambda) ;
130 cfrcheb(deg, deg, coloc, deg, coloc) ;
131 for (int i=0 ; i<n ;i++)
132 res.set(1,i) = coloc[i] ;
133
134 // Limit on the boundaries :
135 val_lim.set(0,0) = pow(r_minus, lambda) ;
136 val_lim.set(0,1) = lambda*pow(r_minus, lambda-1) ;
137 val_lim.set(0,2) = pow(r_plus, lambda) ;
138 val_lim.set(0,3) = lambda*pow(r_plus, lambda-1) ;
139
140 val_lim.set(1,0) = log(r_minus)*pow(r_minus, lambda) ;
141 val_lim.set(1,1) = (1+lambda*log(r_minus))*pow(r_minus, lambda-1) ;
142 val_lim.set(1,2) = log(r_plus)*pow(r_plus, lambda) ;
143 val_lim.set(1,3) = (1+lambda*log(r_plus))*pow(r_plus, lambda-1) ;
144 val_lim /= sqrt(double(2)) ;
145 break ;
146 }
147 case -1:{
148 // Two imaginary solutions :
149 double real_part = (a-b)/2./a ;
150 double imag_part = sqrt(-delta)/2./a ;
151
152 // First SH
153 for (int i=0 ; i<n ; i++)
154 coloc[i] = pow(air[i], real_part)*cos(imag_part*log(air[i])) ;
155 cfrcheb(deg, deg, coloc, deg, coloc) ;
156 for (int i=0 ; i<n ;i++)
157 res.set(0,i) = coloc[i] ;
158
159 // Second SH
160 for (int i=0 ; i<n ; i++)
161 coloc[i] = pow(air[i], real_part)*sin(imag_part*log(air[i])) ;
162 cfrcheb(deg, deg, coloc, deg, coloc) ;
163 for (int i=0 ; i<n ;i++)
164 res.set(1,i) = coloc[i] ;
165
166 // Limit on the boundaries :
167 val_lim.set(0,0) = pow(r_minus, real_part)*cos(imag_part*log(r_minus)) ;
168 val_lim.set(0,1) = (real_part*cos(imag_part*log(r_minus)) -
169 imag_part*sin(imag_part*log(r_minus))) *
170 pow(r_minus, real_part-1) ;
171 val_lim.set(0,2) = pow(r_plus, real_part)*cos(imag_part*log(r_plus)) ;
172 val_lim.set(0,3) = (real_part*cos(imag_part*log(r_plus)) -
173 imag_part*sin(imag_part*log(r_plus))) *
174 pow(r_plus, real_part-1) ;
175
176 val_lim.set(1,0) = pow(r_minus, real_part)*sin(imag_part*log(r_minus)) ;
177 val_lim.set(1,1) = (real_part*sin(imag_part*log(r_minus)) +
178 imag_part*cos(imag_part*log(r_minus))) *
179 pow(r_minus, real_part-1) ;
180 val_lim.set(1,2) = pow(r_plus, real_part)*sin(imag_part*log(r_plus)) ;
181 val_lim.set(1,3) = (real_part*sin(imag_part*log(r_plus)) +
182 imag_part*cos(imag_part*log(r_plus))) *
183 pow(r_plus, real_part-1) ;
184 val_lim /= sqrt(double(2)) ;
185 break ;
186 }
187 default:
188 cout << "What are you doing here ? Get out or I call the police !" << endl;
189 abort() ;
190 break ;
191 }
192
193 delete [] deg ;
194 delete [] coloc ;
195 delete [] air ;
196
197 return res ;
198}
199
200
202
203 // Routines de derivation
205 double, double, double, Tbl&) ;
206 static int nap = 0 ;
207
208 // Premier appel
209 if (nap==0) {
210 nap = 1 ;
211 for (int i=0 ; i<MAX_BASE ; i++) {
212 solh_sec_order_r2[i] = _solh_sec_order_r2_pas_prevu ;
213 }
214 // Les routines existantes
215 solh_sec_order_r2[R_CHEB >> TRA_R] = _solh_sec_order_r2_r_cheb ;
216 }
217
218 Tbl val_lim (2,4) ;
219 val_lim.set_etat_qcq() ;
221 val_lim)) ;
222
223 s_one_minus = val_lim(0,0) ;
224 ds_one_minus = val_lim(0,1) ;
225 s_one_plus = val_lim(0,2) ;
226 ds_one_plus = val_lim(0,3) ;
227
228 s_two_minus = val_lim(1,0) ;
229 ds_two_minus = val_lim(1,1) ;
230 s_two_plus = val_lim(1,2) ;
231 ds_two_plus = val_lim(1,3) ;
232
233 return res ;
234}
235}
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 b .
virtual Tbl get_solh() const
Computes the homogeneous solutions(s).
double a_param
The parameter a .
double c_param
The parameter c .
Basic array class.
Definition tbl.h:161
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:69
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
#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