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 
23 char 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 
72 namespace Lorene {
73 Matrice _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 
88 Matrice _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 
143 Matrice _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  //----------------------------
209 Matrice 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.
Definition: type_parite.h:144
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Lorene prototypes.
Definition: app_hor.h:64