LORENE
helmholtz_plus_mat.C
1 /*
2  * Copyright (c) 1999-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 helmholtz_plus_mat_C[] = "$$" ;
24 
25 /*
26  * $Id: helmholtz_plus_mat.C,v 1.6 2014/10/13 08:53:28 j_novak Exp $
27  * $Log: helmholtz_plus_mat.C,v $
28  * Revision 1.6 2014/10/13 08:53:28 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.5 2014/10/06 15:16:08 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.4 2007/05/06 10:48:12 p_grandclement
35  * Modification of a few operators for the vorton project
36  *
37  * Revision 1.3 2004/01/15 09:15:37 p_grandclement
38  * Modification and addition of the Helmholtz operators
39  *
40  * Revision 1.2 2003/12/11 15:57:26 p_grandclement
41  * include stdlib.h encore ...
42  *
43  * Revision 1.1 2003/12/11 14:48:49 p_grandclement
44  * Addition of ALL (and that is a lot !) the files needed for the general elliptic solver ... UNDER DEVELOPEMENT...
45  *
46  *
47  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/helmholtz_plus_mat.C,v 1.6 2014/10/13 08:53:28 j_novak Exp $
48  *
49  */
50 #include <cstdlib>
51 
52 #include "matrice.h"
53 #include "type_parite.h"
54 #include "proto.h"
55 
56  //-----------------------------------
57  // Routine pour les cas non prevus --
58  //-----------------------------------
59 
60 namespace Lorene {
61 Matrice _helmholtz_plus_mat_pas_prevu(int, int, double, double, double) {
62  cout << "Helmholtz plus : base not implemented..." << endl ;
63  abort() ;
64  exit(-1) ;
65  Matrice res(1, 1) ;
66  return res;
67 }
68 
69 
70 
71  //-------------------------
72  //-- CAS R_CHEBP -----
73  //--------------------------
74 
75 
76 Matrice _helmholtz_plus_mat_r_chebp (int n, int lq, double alpha, double,
77  double masse) {
78  assert (masse > 0) ;
79 
80  Matrice dd(n, n) ;
81  dd.set_etat_qcq() ;
82  Matrice xd(n, n) ;
83  xd.set_etat_qcq() ;
84  Matrice xx(n, n) ;
85  xx.set_etat_qcq() ;
86  Matrice sx2(n, n) ;
87  sx2.set_etat_qcq() ;
88 
89  double* vect = new double[n] ;
90 
91  for (int i=0 ; i<n ; i++) {
92  for (int j=0 ; j<n ; j++)
93  vect[j] = 0 ;
94  vect[i] = 1 ;
95  d2sdx2_1d (n, &vect, R_CHEBP) ;
96  for (int j=0 ; j<n ; j++)
97  dd.set(j, i) = vect[j] ;
98  }
99 
100  for (int i=0 ; i<n ; i++) {
101  for (int j=0 ; j<n ; j++)
102  vect[j] = 0 ;
103  vect[i] = 1 ;
104  sxdsdx_1d (n, &vect, R_CHEBP) ;
105  for (int j=0 ; j<n ; j++)
106  xd.set(j, i) = vect[j] ;
107  }
108 
109  for (int i=0 ; i<n ; i++) {
110  for (int j=0 ; j<n ; j++)
111  vect[j] = 0 ;
112  vect[i] = 1 ;
113  sx2_1d (n, &vect, R_CHEBP) ;
114  for (int j=0 ; j<n ; j++)
115  sx2.set(j, i) = vect[j] ;
116  }
117 
118  for (int i=0 ; i<n ; i++) {
119  for (int j=0 ; j<n ; j++)
120  xx.set(i,j) = 0 ;
121  xx.set(i, i) = 1 ;
122  }
123 
124  delete [] vect ;
125 
126  Matrice res(n, n) ;
127  res = dd+2*xd-lq*(lq+1)*sx2+masse*masse*alpha*alpha*xx ;
128 
129  return res ;
130 }
131 
132 
133  //-------------------------
134  //-- CAS R_CHEB -----
135  //------------------------
136 
137 Matrice _helmholtz_plus_mat_r_cheb (int n, int lq, double alpha, double beta,
138  double masse) {
139 
140 
141 
142  assert (masse > 0) ;
143 
144  double echelle = beta / alpha ;
145 
146  Matrice dd(n, n) ;
147  dd.set_etat_qcq() ;
148  Matrice xd(n, n) ;
149  xd.set_etat_qcq() ;
150  Matrice xx(n, n) ;
151  xx.set_etat_qcq() ;
152 
153  double* vect = new double[n] ;
154 
155  for (int i=0 ; i<n ; i++) {
156  for (int j=0 ; j<n ; j++)
157  vect[j] = 0 ;
158  vect[i] = 1 ;
159  d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
160  vect[i] += masse*masse*alpha*alpha ;
161  for (int j=0 ; j<n ; j++)
162  dd.set(j, i) = vect[j]*echelle*echelle ;
163  }
164 
165  for (int i=0 ; i<n ; i++) {
166  for (int j=0 ; j<n ; j++)
167  vect[j] = 0 ;
168  vect[i] = 1 ;
169  d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
170  vect[i] += masse*masse*alpha*alpha ;
171  multx_1d (n, &vect, R_CHEB) ;
172  for (int j=0 ; j<n ; j++)
173  dd.set(j, i) += 2*echelle*vect[j] ;
174  }
175 
176  for (int i=0 ; i<n ; i++) {
177  for (int j=0 ; j<n ; j++)
178  vect[j] = 0 ;
179  vect[i] = 1 ;
180  d2sdx2_1d (n, &vect, R_CHEB) ; // appel dans le cas fin
181  vect[i] += masse*masse*alpha*alpha ;
182  multx_1d (n, &vect, R_CHEB) ;
183  multx_1d (n, &vect, R_CHEB) ;
184  for (int j=0 ; j<n ; j++)
185  dd.set(j, i) += vect[j] ;
186  }
187 
188  for (int i=0 ; i<n ; i++) {
189  for (int j=0 ; j<n ; j++)
190  vect[j] = 0 ;
191  vect[i] = 1 ;
192  sxdsdx_1d (n, &vect, R_CHEB) ;
193  for (int j=0 ; j<n ; j++)
194  xd.set(j, i) = vect[j]*echelle ;
195  }
196 
197  for (int i=0 ; i<n ; i++) {
198  for (int j=0 ; j<n ; j++)
199  vect[j] = 0 ;
200  vect[i] = 1 ;
201  sxdsdx_1d (n, &vect, R_CHEB) ;
202  multx_1d (n, &vect, R_CHEB) ;
203  for (int j=0 ; j<n ; j++)
204  xd.set(j, i) += vect[j] ;
205  }
206 
207  for (int i=0 ; i<n ; i++) {
208  for (int j=0 ; j<n ; j++)
209  vect[j] = 0 ;
210  vect[i] = 1 ;
211  sx2_1d (n, &vect, R_CHEB) ;
212  for (int j=0 ; j<n ; j++)
213  xx.set(j, i) = vect[j] ;
214  }
215 
216  delete [] vect ;
217 
218  Matrice res(n, n) ;
219  res = dd+2*xd - lq*(lq+1)*xx;
220 
221  return res ;
222 }
223 
224 
225  //--------------------------
226  //- La routine a appeler ---
227  //----------------------------
228 
229 Matrice helmholtz_plus_mat(int n, int lq, double alpha, double beta, double masse,
230  int base_r)
231 {
232 
233  // Routines de derivation
234  static Matrice (*helmholtz_plus_mat[MAX_BASE])(int, int, double, double, double);
235  static int nap = 0 ;
236 
237  // Premier appel
238  if (nap==0) {
239  nap = 1 ;
240  for (int i=0 ; i<MAX_BASE ; i++) {
241  helmholtz_plus_mat[i] = _helmholtz_plus_mat_pas_prevu ;
242  }
243  // Les routines existantes
244  helmholtz_plus_mat[R_CHEB >> TRA_R] = _helmholtz_plus_mat_r_cheb ;
245  helmholtz_plus_mat[R_CHEBP >> TRA_R] = _helmholtz_plus_mat_r_chebp ;
246  }
247 
248  Matrice res(helmholtz_plus_mat[base_r](n, lq, alpha, beta, masse)) ;
249  return res ;
250 }
251 
252 }
Lorene prototypes.
Definition: app_hor.h:64
#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
#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