LORENE
ope_poisson_pseudo_1d_mat.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_pseudo_1d_mat.C,v 1.4 2016/12/05 16:18:13 j_novak Exp $
25  * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_poisson_pseudo_1d/ope_poisson_pseudo_1d_mat.C,v 1.4 2016/12/05 16:18:13 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  //-----------------------------------
37 
38 namespace Lorene {
39 Matrice _poisson_pseudo_1d_mat_pas_prevu(int, int, double, double) {
40  cout << "Operator pas prevu..." << endl ;
41  abort() ;
42  exit(-1) ;
43  Matrice res(1, 1) ;
44  return res;
45 }
46 
47 
48  //-------------------------
49  //-- CAS R_CHEBP -----
50  //--------------------------
51 
52 
53 Matrice _poisson_pseudo_1d_mat_r_chebp (int n, int l, double, double) {
54 
55  Matrice dd(n, n) ;
56  dd.set_etat_qcq() ;
57  Matrice xx(n, n) ;
58  xx.set_etat_qcq() ;
59 
60  double* vect = new double[n] ;
61 
62  for (int i=0 ; i<n ; i++) {
63  for (int j=0 ; j<n ; j++)
64  vect[j] = 0 ;
65  vect[i] = 1 ;
66  d2sdx2_1d (n, &vect, R_CHEBP) ;
67 
68  for (int j=0 ; j<n ; j++)
69  dd.set(j, i) = vect[j] ;
70  }
71 
72  for (int i=0 ; i<n ; i++) {
73  for (int j=0 ; j<n ; j++)
74  vect[j] = 0 ;
75  vect[i] = 1 ;
76  sx2_1d (n, &vect, R_CHEBP) ;
77  for (int j=0 ; j<n ; j++)
78  xx.set(j, i) = vect[j] ;
79  }
80 
81  delete [] vect ;
82 
83  Matrice res(n, n) ;
84  res = dd-l*(l-1)*xx ;
85 
86  return res ;
87 }
88 
89 
90  //------------------------
91  //-- CAS R_CHEBI ----
92  //------------------------
93 
94 
95 Matrice _poisson_pseudo_1d_mat_r_chebi (int n, int l,
96  double, double) {
97 
98  Matrice dd(n, n) ;
99  dd.set_etat_qcq() ;
100  Matrice xx(n, n) ;
101  xx.set_etat_qcq() ;
102 
103  double* vect = new double[n] ;
104 
105  for (int i=0 ; i<n ; i++) {
106  for (int j=0 ; j<n ; j++)
107  vect[j] = 0 ;
108  vect[i] = 1 ;
109  d2sdx2_1d (n, &vect, R_CHEBI) ; // appel dans le cas impair
110  for (int j=0 ; j<n ; j++)
111  dd.set(j, i) = vect[j] ;
112  }
113 
114  for (int i=0 ; i<n ; i++) {
115  for (int j=0 ; j<n ; j++)
116  vect[j] = 0 ;
117  vect[i] = 1 ;
118  sx2_1d (n, &vect, R_CHEBI) ;
119  for (int j=0 ; j<n ; j++)
120  xx.set(j, i) = vect[j] ;
121  }
122 
123  delete [] vect ;
124 
125  Matrice res(n, n) ;
126  res = dd-l*(l-1)*xx ;
127  return res ;
128 }
129 
130  //-------------------------
131  //-- CAS R_CHEB -----
132  //-----------------------
133 
134 
135 Matrice _poisson_pseudo_1d_mat_r_cheb (int n, int l,
136  double alf, double bet) {
137 
138  double echelle = bet / alf ;
139 
140 
141  Matrice dd(n, n) ;
142  dd.set_etat_qcq() ;
143  Matrice xx(n, n) ;
144  xx.set_etat_qcq() ;
145 
146  double* vect = new double[n] ;
147 
148  for (int i=0 ; i<n ; i++) {
149  for (int j=0 ; j<n ; j++)
150  vect[j] = 0 ;
151  vect[i] = 1 ;
152  d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
153  for (int j=0 ; j<n ; j++)
154  dd.set(j, i) = vect[j]*echelle*echelle ;
155  }
156 
157  for (int i=0 ; i<n ; i++) {
158  for (int j=0 ; j<n ; j++)
159  vect[j] = 0 ;
160  vect[i] = 1 ;
161  d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
162  multx_1d (n, &vect, R_CHEB) ;
163  for (int j=0 ; j<(n>i+1 ? i+1 : n) ; j++)
164  dd.set(j, i) += 2*echelle*vect[j] ;
165  }
166 
167  for (int i=0 ; i<n ; i++) {
168  for (int j=0 ; j<n ; j++)
169  vect[j] = 0 ;
170  vect[i] = 1 ;
171  d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
172  multx_1d (n, &vect, R_CHEB) ;
173  multx_1d (n, &vect, R_CHEB) ;
174  for (int j=0 ; j<(n>i+1 ? i+1 : n) ; j++)
175  dd.set(j, i) += vect[j] ;
176  }
177 
178  for (int i=0 ; i<n ; i++) {
179  for (int j=0 ; j<n ; j++)
180  vect[j] = 0 ;
181  vect[i] = 1 ;
182  sx2_1d (n, &vect, R_CHEB) ;
183  for (int j=0 ; j<n ; j++)
184  xx.set(j, i) = vect[j] ;
185  }
186 
187  delete [] vect ;
188 
189  Matrice res(n, n) ;
190  res = dd-l*(l-1)*xx ;
191  return res ;
192 }
193 
194 
196  if (ope_mat != 0x0)
197  delete ope_mat ;
198 
199  // Routines de derivation
200  static Matrice (*poisson_pseudo_1d_mat[MAX_BASE])(int, int, double, double);
201  static int nap = 0 ;
202 
203  // Premier appel
204  if (nap==0) {
205  nap = 1 ;
206  for (int i=0 ; i<MAX_BASE ; i++) {
207  poisson_pseudo_1d_mat[i] = _poisson_pseudo_1d_mat_pas_prevu ;
208  }
209  // Les routines existantes
210  poisson_pseudo_1d_mat[R_CHEB >> TRA_R] = _poisson_pseudo_1d_mat_r_cheb ;
211  poisson_pseudo_1d_mat[R_CHEBP >> TRA_R] = _poisson_pseudo_1d_mat_r_chebp ;
212  poisson_pseudo_1d_mat[R_CHEBI >> TRA_R] = _poisson_pseudo_1d_mat_r_chebi ;
213  }
214  ope_mat = new Matrice(poisson_pseudo_1d_mat[base_r](nr, l_quant,
215  alpha, beta)) ;
216 }
217 }
double alpha
Parameter of the associated mapping.
double beta
Parameter of the associated mapping.
Lorene prototypes.
Definition: app_hor.h:67
Matrice * ope_mat
Pointer on the matrix representation of the operator.
int base_r
Radial basis of decomposition.
#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
virtual void do_ope_mat() const
Computes the matrix of the operator.
Matrix handling.
Definition: matrice.h:152
int nr
Number of radial points.
#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