LORENE
ope_vorton_solp.C
1 /*
2  * Copyright (c) 2003 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_vorton_solp.C,v 1.5 2016/12/05 16:18:14 j_novak Exp $
25  * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_vorton/ope_vorton_solp.C,v 1.5 2016/12/05 16:18:14 j_novak Exp $
26  *
27  */
28 #include <cmath>
29 #include <cstdlib>
30 
31 #include "proto.h"
32 #include "ope_elementary.h"
33 
34 
35  //------------------------------------
36  // Cl version Tbl -> Tbl --
37  //------------------------------------
38 namespace Lorene {
39 Tbl _cl_vorton_pas_prevu (const Tbl &so, int) {
40 
41  cout << "Linear combination for vorton not implemented..." << endl ;
42  abort() ;
43  exit(-1) ;
44  return so;
45 }
46 
47  //-------------------
48  //-- R_CHEB -------
49  //--------------------
50 Tbl _cl_vorton_r_cheb (const Tbl& source, int) {
51 
52  int n = source.get_dim(0) ;
53 
54  Tbl barre(source) ;
55  int dirac = 1 ;
56  for (int i=0 ; i<n-2 ; i++) {
57  barre.set(i) = ((1+dirac)*source(i)-source(i+2))
58  /(i+1) ;
59  if (i==0) dirac = 0 ;
60  }
61 
62  Tbl res(barre) ;
63  for (int i=0 ; i<n-4 ; i++)
64  res.set(i) = barre(i)-barre(i+2) ;
65 
66  return res ;
67 }
68 
69  //-------------------
70  //-- R_CHEBU -----
71  //-------------------
72 
73 Tbl _cl_vorton_r_chebu_trois (const Tbl &source) {
74  int n = source.get_dim(0) ;
75  Tbl barre(source) ;
76  int dirac = 1 ;
77  for (int i=0 ; i<n-2 ; i++) {
78  barre.set(i) = ((1+dirac)*source(i)-source(i+2)) ;
79  if (i==0) dirac = 0 ;
80  }
81 
82  Tbl tilde(barre) ;
83  for (int i=0 ; i<n-4 ; i++)
84  tilde.set(i) = (barre(i)-barre(i+2)) ;
85 
86  Tbl res(tilde) ;
87  for (int i=0 ; i<n-4 ; i++)
88  res.set(i) = (tilde(i)+tilde(i+1)) ;
89  return res ;
90 }
91 
92 Tbl _cl_vorton_r_chebu (const Tbl& source, int puis) {
93  int n = source.get_dim(0) ;
94  Tbl res (n) ;
95  res.set_etat_qcq() ;
96 
97  switch (puis) {
98  case 3 :
99  res = _cl_vorton_r_chebu_trois(source) ;
100  break ;
101  default :
102  abort() ;
103  exit(-1) ;
104  }
105 return res ;
106 }
107 
108  //----------------------------
109  //- Routine a appeler ---
110  //------------------------------
111 
112 Tbl cl_vorton (const Tbl &source, int puis, int base_r) {
113 
114  // Routines de derivation
115  static Tbl (*cl_vorton[MAX_BASE])(const Tbl &, int) ;
116  static int nap = 0 ;
117 
118  // Premier appel
119  if (nap==0) {
120  nap = 1 ;
121  for (int i=0 ; i<MAX_BASE ; i++) {
122  cl_vorton[i] = _cl_vorton_pas_prevu ;
123  }
124  // Les routines existantes
125  cl_vorton[R_CHEB >> TRA_R] = _cl_vorton_r_cheb ;
126  cl_vorton[R_CHEBU >> TRA_R] = _cl_vorton_r_chebu ;
127  }
128 
129  Tbl res(cl_vorton[base_r](source, puis)) ;
130  return res ;
131 }
132 
133 
134  //*******************************
135  // CALCUL SP proprement parler
136  //*******************************
137 
138  //------------------------------------
139  // Routine pour les cas non prevus --
140  //------------------------------------
141 Tbl _solp_vorton_pas_prevu (const Matrice &, const Matrice &,
142  const Tbl &, double, double, int) {
143  cout << " Solution particuliere pas prevue in sec_order..... : "<< endl ;
144  abort() ;
145  exit(-1) ;
146  Tbl res(1) ;
147  return res;
148 }
149  //-------------------
150  //-- R_CHEBU -----
151  //-------------------
152 
153 Tbl _solp_vorton_r_chebu_trois (const Matrice &lap, const Matrice &nondege,
154  const Tbl &source, double alpha) {
155 
156  int n = lap.get_dim(0) ;
157  int dege = n-nondege.get_dim(0) ;;
158  assert ((dege==2) || (dege==1)) ;
159 
160  Tbl source_aux (alpha*cl_vorton (source, 3, R_CHEBU)) ;
161 
162  Tbl so(n-dege) ;
163  so.set_etat_qcq() ;
164  for (int i=0 ; i<n-dege ; i++)
165  so.set(i) = source_aux(i) ;
166 
167  Tbl auxi(nondege.inverse(so)) ;
168 
169  Tbl res(n) ;
170  res.set_etat_qcq() ;
171  for (int i=dege ; i<n ; i++)
172  res.set(i) = auxi(i-dege) ;
173 
174  for (int i=0 ; i<dege ; i++)
175  res.set(i) = 0 ;
176 
177  if (dege==2) {
178  double somme = 0 ;
179  for (int i=0 ; i<n ; i++)
180  somme -= res(i) ;
181  res.set(0) = somme ;
182  }
183 
184  return res ;
185 }
186 
187 
188 Tbl _solp_vorton_r_chebu (const Matrice &lap, const Matrice &nondege,
189  const Tbl &source,double alpha, double, int puis) {
190  int n = source.get_dim(0) ;
191  Tbl res (n) ;
192  res.set_etat_qcq() ;
193 
194  switch (puis) {
195  case 3 :
196  res = _solp_vorton_r_chebu_trois(lap, nondege, source, alpha) ;
197  break ;
198  default :
199  abort() ;
200  exit(-1) ;
201  }
202 return res ;
203 }
204 
205  //-------------------
206  //-- R_CHEB -----
207  //-------------------
208 
209 Tbl _solp_vorton_r_cheb (const Matrice &lap, const Matrice &nondege,
210  const Tbl &source,double alpha, double beta, int dz) {
211 
212  int n = lap.get_dim(0) ;
213  int dege = n-nondege.get_dim(0) ;
214  assert (dege ==2) ;
215 
216  Tbl source_aux(source*alpha*alpha) ;
217  Tbl xso(source_aux) ;
218  Tbl xxso(source_aux) ;
219  multx_1d(n, &xso.t, R_CHEB) ;
220  multx_1d(n, &xxso.t, R_CHEB) ;
221  multx_1d(n, &xxso.t, R_CHEB) ;
222  source_aux = beta*beta/alpha/alpha*source_aux+2*beta/alpha*xso+xxso ;
223  source_aux = cl_vorton (source_aux, dz, R_CHEB) ;
224 
225  Tbl so(n-dege) ;
226  so.set_etat_qcq() ;
227  for (int i=0 ; i<n-dege ; i++)
228  so.set(i) = source_aux(i) ;
229 
230  Tbl auxi(nondege.inverse(so)) ;
231 
232  Tbl res(n) ;
233  res.set_etat_qcq() ;
234  for (int i=dege ; i<n ; i++)
235  res.set(i) = auxi(i-dege) ;
236 
237  for (int i=0 ; i<dege ; i++)
238  res.set(i) = 0 ;
239  return res ;
240 }
241 
242 
243 
244 Tbl Ope_vorton::get_solp (const Tbl& so) const {
245 
246  if (non_dege == 0x0)
247  do_non_dege() ;
248 
249  // Routines de derivation
250  static Tbl (*solp_vorton[MAX_BASE]) (const Matrice&, const Matrice&,
251  const Tbl&, double, double, int) ;
252  static int nap = 0 ;
253 
254  // Premier appel
255  if (nap==0) {
256  nap = 1 ;
257  for (int i=0 ; i<MAX_BASE ; i++) {
258  solp_vorton[i] = _solp_vorton_pas_prevu ;
259  }
260  // Les routines existantes
261  solp_vorton[R_CHEB >> TRA_R] = _solp_vorton_r_cheb ;
262  solp_vorton[R_CHEBU >> TRA_R] = _solp_vorton_r_chebu ;
263  }
264 
265  Tbl res(solp_vorton[base_r] (*ope_mat, *non_dege, so, alpha, beta, dzpuis)) ;
266  Tbl valeurs (val_solp (res, alpha, base_r)) ;
267 
268  sp_plus = valeurs(0) ;
269  sp_minus = valeurs(1) ;
270  dsp_plus = valeurs(2) ;
271  dsp_minus = valeurs(3) ;
272 
273  return res ;
274 }
275 }
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.
double dsp_minus
Value of the derivative of the particular solution at the inner boundary.
int dzpuis
the associated dzpuis, if in the compactified domain.
double sp_minus
Value of the particular solution at the inner boundary.
int base_r
Radial basis of decomposition.
double sp_plus
Value of the particular solution at the outer boundary.
virtual void do_non_dege() const
Computes the non-degenerated matrix of the operator.
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
double dsp_plus
Value of the derivative of the particular solution at the outer boundary.
Matrix handling.
Definition: matrice.h:152
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:423
virtual Tbl get_solp(const Tbl &so) const
Computes the particular solution, given the source so .
int get_dim(int i) const
Returns the dimension of the matrix.
Definition: matrice.C:263
Basic array class.
Definition: tbl.h:164
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
Matrice * non_dege
Pointer on the non-degenerated matrix of the operator.
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166