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
21char ope_poisson_pseudo_1d_mat_C[] = "$Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_poisson_pseudo_1d/ope_poisson_pseudo_1d_mat.C,v 1.3 2014/10/13 08:53:35 j_novak Exp $" ;
22
23/*
24 * $Id: ope_poisson_pseudo_1d_mat.C,v 1.3 2014/10/13 08:53:35 j_novak Exp $
25 * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_poisson_pseudo_1d/ope_poisson_pseudo_1d_mat.C,v 1.3 2014/10/13 08:53:35 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
38namespace Lorene {
39Matrice _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
53Matrice _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
95Matrice _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
135Matrice _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
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 }
215 alpha, beta)) ;
216}
217}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Matrix handling.
Definition matrice.h:152
Matrice * ope_mat
Pointer on the matrix representation of the operator.
double beta
Parameter of the associated mapping.
double alpha
Parameter of the associated mapping.
int base_r
Radial basis of decomposition.
int nr
Number of radial points.
virtual void do_ope_mat() const
Computes the matrix of the operator.
#define MAX_BASE
Nombre max. de bases differentes.
#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