LORENE
xdsdx_mat.C
1/*
2 * Copyright (c) 2000-2001 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 as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23char xdsdx_mat_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/xdsdx_mat.C,v 1.5 2014/10/13 08:53:31 j_novak Exp $" ;
24
25/*
26 * $Id: xdsdx_mat.C,v 1.5 2014/10/13 08:53:31 j_novak Exp $
27 * $Log: xdsdx_mat.C,v $
28 * Revision 1.5 2014/10/13 08:53:31 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.4 2014/10/06 15:16:11 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.3 2007/06/21 20:07:16 k_taniguchi
35 * nmax increased to 200
36 *
37 * Revision 1.2 2002/10/16 14:37:12 j_novak
38 * Reorganization of #include instructions of standard C++, in order to
39 * use experimental version 3 of gcc.
40 *
41 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
42 * LORENE
43 *
44 * Revision 2.0 2000/03/16 16:23:17 phil
45 * *** empty log message ***
46 *
47 *
48 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/xdsdx_mat.C,v 1.5 2014/10/13 08:53:31 j_novak Exp $
49 *
50 */
51
52/*
53 * Routine renvoyan la matrice de l'operateur x f' = s
54 * Pour l != 1 le resultat est donne en s est donne en Chebyshev et
55 * f en Gelerkin (T_i + T_{i+1} pour l pair et (2*i+3)T_i + (2*i+1)T_{i+1} pour
56 * l impair.
57 * Pour l=1 pas de probleme de singularite on reste donc en Chebyshev.
58 */
59
60//fichiers includes
61#include <cstdio>
62#include <cstdlib>
63
64#include "matrice.h"
65#include "type_parite.h"
66#include "proto.h"
67
68 //-----------------------------------
69 // Routine pour les cas non prevus --
70 //-----------------------------------
71
72namespace Lorene {
73Matrice _xdsdx_mat_pas_prevu(int n, int) {
74 cout << "xdsdx_mat pas prevu..." << endl ;
75 cout << "n : " << n << endl ;
76 abort() ;
77 exit(-1) ;
78 Matrice res(1, 1) ; // On ne passe jamais ici de toute facon !
79 return res;
80}
81
82
83 //-------------------------
84 //-- CAS R_CHEBP -----
85 //--------------------------
86
87
88Matrice _xdsdx_mat_r_chebp (int n, int) {
89 const int nmax = 200 ;// Nombre de Matrices stockees
90 static Matrice* tab[nmax] ; // les matrices calculees
91 static int nb_dejafait = 0 ; // nbre de matrices calculees
92 static int nr_dejafait[nmax] ;
93
94 int indice = -1 ;
95
96 // On determine si la matrice a deja ete calculee :
97 for (int conte=0 ; conte<nb_dejafait ; conte ++)
98 if (nr_dejafait[conte] == n)
99 indice = conte ;
100
101 // Calcul a faire :
102 if (indice == -1) {
103 if (nb_dejafait >= nmax) {
104 cout << "_laplacien_nul_mat_r_chebp : trop de matrices" << endl ;
105 abort() ;
106 exit (-1) ;
107 }
108
109 nr_dejafait[nb_dejafait] = n ;
110
111 Matrice res(n-1, n-1) ;
112 res.set_etat_qcq() ;
113
114 double* xdsdx = new double[n] ;
115
116 for (int i=0 ; i<n-1 ; i++) {
117 for (int j=0 ; j<n ; j++)
118 xdsdx[j] = 0 ;
119 xdsdx[i] = 1 ;
120 xdsdx[i+1] = 1 ;
121 xdsdx_1d (n, &xdsdx, R_CHEBP) ;
122
123 for (int j=0 ; j<n-1 ; j++)
124 res.set(j, i) = xdsdx[j] ;
125 }
126 delete [] xdsdx ;
127 tab[nb_dejafait] = new Matrice(res) ;
128 nb_dejafait ++ ;
129 return res ;
130 }
131
132 else
133 return *tab[indice] ;
134}
135
136
137
138 //------------------------
139 //-- CAS R_CHEBI ----
140 //------------------------
141
142
143Matrice _xdsdx_mat_r_chebi (int n, int l) {
144 const int nmax = 200 ;// Nombre de Matrices stockees
145 static Matrice* tab[nmax] ; // les matrices calculees
146 static int nb_dejafait = 0 ; // nbre de matrices calculees
147 static int nr_dejafait[nmax] ;
148 static int nl_dejafait[nmax] ;
149
150 int indice = -1 ;
151 // On separe les cas l=1 et l != 1
152 int indic_l = (l == 1) ? 1 : 2 ;
153
154 // On determine si la matrice a deja ete calculee :
155 for (int conte=0 ; conte<nb_dejafait ; conte ++)
156 if ((nr_dejafait[conte] == n) && (nl_dejafait[conte] == indic_l))
157 indice = conte ;
158
159 // Calcul a faire :
160 if (indice == -1) {
161 if (nb_dejafait >= nmax) {
162 cout << "_laplacien_nul_mat_r_chebp : trop de matrices" << endl ;
163 abort() ;
164 exit (-1) ;
165 }
166
167 nr_dejafait[nb_dejafait] = n ;
168 nl_dejafait[nb_dejafait] = indic_l ;
169
170 // non degenere pour l=1
171 int taille = (l==1) ? n : n-1 ;
172 Matrice res(taille, taille) ;
173 res.set_etat_qcq() ;
174
175 double* xdsdx = new double[n] ;
176
177 for (int i=0 ; i<taille ; i++) {
178 for (int j=0 ; j<n ; j++)
179 xdsdx[j] = 0 ;
180
181 // Gelerkin ou Chebyshev ?
182 if (taille == n) {
183 xdsdx[i] = 1 ;
184 }
185 else {
186 xdsdx[i] = 2*i+3 ;
187 xdsdx[i+1] = 2*i+1 ;
188 }
189 xdsdx_1d (n, &xdsdx, R_CHEBI) ; // appel dans le cas impair
190
191 for (int j=0 ; j<taille ; j++)
192 res.set(j, i) = xdsdx[j] ;
193 }
194
195 delete [] xdsdx ;
196 tab[nb_dejafait] = new Matrice(res) ;
197 nb_dejafait ++ ;
198 return res ;
199 }
200
201 else
202 return *tab[indice] ;
203}
204
205
206 //--------------------------
207 //- La routine a appeler ---
208 //----------------------------
209Matrice xdsdx_mat(int n, int l, int base_r)
210{
211
212 // Routines de derivation
213 static Matrice (*xdsdx_mat[MAX_BASE])(int, int) ;
214 static int nap = 0 ;
215
216 // Premier appel
217 if (nap==0) {
218 nap = 1 ;
219 for (int i=0 ; i<MAX_BASE ; i++) {
220 xdsdx_mat[i] = _xdsdx_mat_pas_prevu ;
221 }
222 // Les routines existantes
223 xdsdx_mat[R_CHEBP >> TRA_R] = _xdsdx_mat_r_chebp ;
224 xdsdx_mat[R_CHEBI >> TRA_R] = _xdsdx_mat_r_chebi ;
225 }
226
227 Matrice res(xdsdx_mat[base_r](n, l)) ;
228 return res ;
229}
230
231}
#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_CHEBP
base de Cheb. paire (rare) seulement
Lorene prototypes.
Definition app_hor.h:64