LORENE
ope_poisson_2d_solh.C
1 /*
2  * Copyright (c) 2004 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 
21 
22 
23 /*
24  * $Id: ope_poisson_2d_solh.C,v 1.4 2016/12/05 16:18:12 j_novak Exp $
25  * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_poisson_2d/ope_poisson_2d_solh.C,v 1.4 2016/12/05 16:18:12 j_novak Exp $
26  *
27  */
28 #include <cmath>
29 #include <cstdlib>
30 
31 #include "proto.h"
32 #include "ope_elementary.h"
33 
34 
35  //------------------------------------
36  // Routine pour les cas non prevus --
37  //------------------------------------
38 namespace Lorene {
39 Tbl _solh_poisson_2d_pas_prevu (int, int, double, double, Tbl&) {
40 
41  cout << " Solution homogene pas prevue ..... : "<< endl ;
42  exit(-1) ;
43  Tbl res(1) ;
44  return res;
45 }
46 
47 
48  //-------------------
49  //-- R_CHEB ------
50  //-------------------
51 
52 Tbl _solh_poisson_2d_r_cheb (int n, int l, double alpha, double beta,
53  Tbl& val_lim) {
54 
55 
56  double echelle = beta / alpha ;
57 
58  // CASE l != 0
59  if (l > 0) {
60  val_lim.set(0,0) = pow(echelle-1, double(l)) ;
61  val_lim.set(0,1) = double(l) * pow(echelle-1, double(l-1))/alpha ;
62  val_lim.set(0,2) = pow(echelle+1, double(l)) ;
63  val_lim.set(0,3) = double(l) * pow(echelle+1, double(l-1))/alpha ;
64 
65 
66  val_lim.set(1,0) = pow(echelle-1, -double(l)) ;
67  val_lim.set(1,1) = -double(l) * pow(echelle-1, -double(l+1))/alpha ;
68  val_lim.set(1,2) = pow(echelle+1, -double(l)) ;
69  val_lim.set(1,3) = -double(l) * pow(echelle+1, -double(l+1))/alpha ;
70  }
71  // CASE l =0
72  else {
73  val_lim.set(0,0) = 1. ;
74  val_lim.set(0,1) = 0. ;
75  val_lim.set(0,2) = 1. ;
76  val_lim.set(0,3) = 0. ;
77 
78  val_lim.set(1,0) = log(echelle-1) ;
79  val_lim.set(1,1) = 1./(echelle-1)/alpha ;
80  val_lim.set(1,2) = log(echelle+1) ;
81  val_lim.set(1,3) = 1./(echelle+1)/alpha ;
82  }
83 
84  const int nmax = 200 ; // Nombre de Tbl stockes
85  static Tbl* tab[nmax] ; // les Tbl calcules
86  static int nb_dejafait = 0 ; // nbre de Tbl calcules
87  static int l_dejafait[nmax] ;
88  static int nr_dejafait[nmax] ;
89  static double vieux_echelle = 0;
90 
91  // Si on a change l'echelle : on detruit tout :
92  if (vieux_echelle != echelle) {
93  for (int i=0 ; i<nb_dejafait ; i++) {
94  l_dejafait[i] = -1 ;
95  nr_dejafait[i] = -1 ;
96  delete tab[i] ;
97  }
98  nb_dejafait = 0 ;
99  vieux_echelle = echelle ;
100  }
101 
102  int indice = -1 ;
103 
104  // On determine si la matrice a deja ete calculee :
105  for (int conte=0 ; conte<nb_dejafait ; conte ++)
106  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
107  indice = conte ;
108 
109  // Calcul a faire :
110  if (indice == -1) {
111  if (nb_dejafait >= nmax) {
112  cout << "_solh_poisson_2d_r_cheb : trop de Tbl" << endl ;
113  abort() ;
114  exit (-1) ;
115  }
116 
117 
118  l_dejafait[nb_dejafait] = l ;
119  nr_dejafait[nb_dejafait] = n ;
120 
121  Tbl res(2, n) ;
122  res.set_etat_qcq() ;
123  double* coloc = new double[n] ;
124 
125  int * deg = new int[3] ;
126  deg[0] = 1 ;
127  deg[1] = 1 ;
128  deg[2] = n ;
129 
130  //Construction de la premiere solution homogene :
131  // cad celle polynomiale.
132 
133  for (int i=0 ; i<n ; i++)
134  if (l!=0)
135  coloc[i] = pow(echelle-cos(M_PI*i/(n-1)), double(l)) ;
136  else
137  coloc[i] = 1 ;
138 
139  cfrcheb(deg, deg, coloc, deg, coloc) ;
140  for (int i=0 ; i<n ;i++)
141  res.set(0, i) = coloc[i] ;
142 
143  // construction de la seconde solution homogene :
144  // cad celle fractionnelle (ou log dans le cas l==0)
145  for (int i=0 ; i<n ; i++)
146  if (l != 0)
147  coloc[i] = pow(echelle-cos(M_PI*i/(n-1)), -double(l)) ;
148  else
149  coloc[i] = log(echelle-cos(M_PI*i/(n-1))) ;
150 
151  cfrcheb(deg, deg, coloc, deg, coloc) ;
152  for (int i=0 ; i<n ;i++)
153  res.set(1, i) = coloc[i] ;
154 
155 
156  delete [] coloc ;
157  delete [] deg ;
158  tab[nb_dejafait] = new Tbl(res) ;
159  nb_dejafait ++ ;
160  return res ;
161  }
162 
163  else return *tab[indice] ;
164 }
165 
166  //-------------------
167  //-- R_CHEBP ------
168  //-------------------
169 
170 Tbl _solh_poisson_2d_r_chebp (int n, int l, double alpha,
171  double, Tbl& val_lim) {
172 
173  val_lim.set(0,0) = (l!=0) ? 1 : 0 ;
174  val_lim.set(0,1) = (l!=1) ? 0 : 1 ;
175  val_lim.set(0,2) = 1. ;
176  val_lim.set(0,3) = double(l)/alpha ;
177 
178  const int nmax = 200 ; // Nombre de Tbl stockes
179  static Tbl* tab[nmax] ; // les Tbl calcules
180  static int nb_dejafait = 0 ; // nbre de Tbl calcules
181  static int l_dejafait[nmax] ;
182  static int nr_dejafait[nmax] ;
183 
184  int indice = -1 ;
185 
186  // On determine si la matrice a deja ete calculee :
187  for (int conte=0 ; conte<nb_dejafait ; conte ++)
188  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
189  indice = conte ;
190 
191  // Calcul a faire :
192  if (indice == -1) {
193  if (nb_dejafait >= nmax) {
194  cout << "_solh_poisson_2d_r_chebp : trop de Tbl" << endl ;
195  abort() ;
196  exit (-1) ;
197  }
198 
199 
200  l_dejafait[nb_dejafait] = l ;
201  nr_dejafait[nb_dejafait] = n ;
202 
203  assert (div(l, 2).rem ==0) ;
204 
205  Tbl res(n) ;
206  res.set_etat_qcq() ;
207  double* coloc = new double[n] ;
208 
209  int * deg = new int[3] ;
210  deg[0] = 1 ;
211  deg[1] = 1 ;
212  deg[2] = n ;
213 
214  for (int i=0 ; i<n ; i++)
215  if (l!=0)
216  coloc[i] = pow(sin(M_PI*i/2/(n-1)), double(l)) ;
217  else
218  coloc[i] = 1 ;
219 
220  cfrchebp(deg, deg, coloc, deg, coloc) ;
221  for (int i=0 ; i<n ;i++)
222  res.set(i) = coloc[i] ;
223 
224  delete [] coloc ;
225  delete [] deg ;
226  tab[nb_dejafait] = new Tbl(res) ;
227  nb_dejafait ++ ;
228  return res ;
229  }
230 
231  else return *tab[indice] ;
232 }
233 
234 
235  //-------------------
236  //-- R_CHEBI -----
237  //-------------------
238 
239 Tbl _solh_poisson_2d_r_chebi (int n, int l,
240  double alpha, double, Tbl& val_lim) {
241 
242  val_lim.set(0,0) = 0 ;
243  val_lim.set(0,1) = (l!=1) ? 0 : 1 ;
244  val_lim.set(0,2) = 1. ;
245  val_lim.set(0,3) = double(l)/alpha ;
246 
247  const int nmax = 200 ; // Nombre de Tbl stockes
248  static Tbl* tab[nmax] ; // les Tbl calcules
249  static int nb_dejafait = 0 ; // nbre de Tbl calcules
250  static int l_dejafait[nmax] ;
251  static int nr_dejafait[nmax] ;
252 
253  int indice = -1 ;
254 
255  // On determine si la matrice a deja ete calculee :
256  for (int conte=0 ; conte<nb_dejafait ; conte ++)
257  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
258  indice = conte ;
259 
260  // Calcul a faire :
261  if (indice == -1) {
262  if (nb_dejafait >= nmax) {
263  cout << "_solh_r_chebi : trop de Tbl" << endl ;
264  abort() ;
265  exit (-1) ;
266  }
267 
268 
269  l_dejafait[nb_dejafait] = l ;
270  nr_dejafait[nb_dejafait] = n ;
271 
272 
273  assert (div(l, 2).rem == 1) ;
274 
275  Tbl res(n) ;
276  res.set_etat_qcq() ;
277  double* coloc = new double[n] ;
278 
279  int * deg = new int[3] ;
280  deg[0] = 1 ;
281  deg[1] = 1 ;
282  deg[2] = n ;
283 
284  for (int i=0 ; i<n ; i++)
285  coloc[i] = pow(sin(M_PI*i/2/(n-1)), double(l)) ;
286 
287  cfrchebi(deg, deg, coloc, deg, coloc) ;
288  for (int i=0 ; i<n ;i++)
289  res.set(i) = coloc[i] ;
290 
291  delete [] coloc ;
292  delete [] deg ;
293  tab[nb_dejafait] = new Tbl(res) ;
294  nb_dejafait ++ ;
295  return res ;
296  }
297 
298  else return *tab[indice] ;
299 }
300 
301 
302 
303  //-------------------
304  //-- R_CHEBU -----
305  //-------------------
306 
307 Tbl _solh_poisson_2d_r_chebu (int n, int l, double alpha,
308  double, Tbl& val_lim) {
309 
310 
311  if (l==0) {
312  cout << "Case l=0 in 2D Poisson not defined in the external compactified domain..." << endl ;
313  abort() ;
314  }
315 
316  val_lim.set(0,0) = pow(-2., double(l)) ;
317  val_lim.set(0,1) = -double(l)*alpha*pow(-2, double(l+1.)) ;
318  val_lim.set(0,2) = 0 ;
319  val_lim.set(0,3) = 0 ;
320 
321  const int nmax = 200 ; // Nombre de Tbl stockes
322  static Tbl* tab[nmax] ; // les Tbl calcules
323  static int nb_dejafait = 0 ; // nbre de Tbl calcules
324  static int l_dejafait[nmax] ;
325  static int nr_dejafait[nmax] ;
326 
327  int indice = -1 ;
328 
329  // On determine si la matrice a deja ete calculee :
330  for (int conte=0 ; conte<nb_dejafait ; conte ++)
331  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
332  indice = conte ;
333 
334  // Calcul a faire :
335  if (indice == -1) {
336  if (nb_dejafait >= nmax) {
337  cout << "_solh_poisson_2d_r_chebu : trop de Tbl" << endl ;
338  abort() ;
339  exit (-1) ;
340  }
341 
342 
343  l_dejafait[nb_dejafait] = l ;
344  nr_dejafait[nb_dejafait] = n ;
345 
346 
347  // assert (l < n-1) ;
348  Tbl res(n) ;
349  res.set_etat_qcq() ;
350  double* coloc = new double[n] ;
351 
352  int * deg = new int[3] ;
353  deg[0] = 1 ;
354  deg[1] = 1 ;
355  deg[2] = n ;
356 
357  for (int i=0 ; i<n ; i++)
358  coloc[i] = pow(-1-cos(M_PI*i/(n-1)), double(l)) ;
359 
360  cfrcheb(deg, deg, coloc, deg, coloc) ;
361  for (int i=0 ; i<n ;i++)
362  res.set(i) = coloc[i] ;
363 
364  delete [] coloc ;
365  delete [] deg ;
366  tab[nb_dejafait] = new Tbl(res) ;
367  nb_dejafait ++ ;
368  return res ;
369  }
370 
371  else return *tab[indice] ;
372 }
373 
374 
376 
377  // Routines de derivation
378  static Tbl (*solh_poisson_2d[MAX_BASE]) (int, int, double, double, Tbl&) ;
379  static int nap = 0 ;
380 
381  // Premier appel
382  if (nap==0) {
383  nap = 1 ;
384  for (int i=0 ; i<MAX_BASE ; i++) {
385  solh_poisson_2d[i] = _solh_poisson_2d_pas_prevu ;
386  }
387  // Les routines existantes
388  solh_poisson_2d[R_CHEB >> TRA_R] = _solh_poisson_2d_r_cheb ;
389  solh_poisson_2d[R_CHEBP >> TRA_R] = _solh_poisson_2d_r_chebp ;
390  solh_poisson_2d[R_CHEBI >> TRA_R] = _solh_poisson_2d_r_chebi ;
391  solh_poisson_2d[R_CHEBU >> TRA_R] = _solh_poisson_2d_r_chebu ;
392  }
393 
394  Tbl val_lim (2,4) ;
395  val_lim.set_etat_qcq() ;
396  Tbl res(solh_poisson_2d[base_r](nr,l_quant, alpha, beta, val_lim)) ;
397 
398  s_one_minus = val_lim(0,0) ;
399  ds_one_minus = val_lim(0,1) ;
400  s_one_plus = val_lim(0,2) ;
401  ds_one_plus = val_lim(0,3) ;
402 
403  s_two_minus = val_lim(1,0) ;
404  ds_two_minus = val_lim(1,1) ;
405  s_two_plus = val_lim(1,2) ;
406  ds_two_plus = val_lim(1,3) ;
407 
408  return res ;
409 }
410 }
double alpha
Parameter of the associated mapping.
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
double s_one_minus
Value of the first homogeneous solution at the inner boundary.
double ds_two_minus
Value of the derivative of the second homogeneous solution at the inner boundary. ...
double beta
Parameter of the associated mapping.
Lorene prototypes.
Definition: app_hor.h:67
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
double ds_two_plus
Value of the derivative of the second homogeneous solution at the outer boundary. ...
int l_quant
quantum number
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:97
int base_r
Radial basis of decomposition.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
double ds_one_plus
Value of the derivative of the first homogeneous solution at the outer boundary.
#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
virtual Tbl get_solh() const
Computes the homogeneous solutions(s).
double s_two_minus
Value of the second homogeneous solution at the inner boundary.
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
double s_one_plus
Value of the first homogeneous solution at the outer boundary.
int nr
Number of radial points.
double ds_one_minus
Value of the derivative of the first homogeneous solution at the inner boundary.
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:72
Basic array class.
Definition: tbl.h:164
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
double s_two_plus
Value of the second homogeneous solution at the outer boundary.
#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