LORENE
ope_poisson_pseudo_1d_solp.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 
21 
22 
23 /*
24  * $Id: ope_poisson_pseudo_1d_solp.C,v 1.4 2016/12/05 16:18:13 j_novak Exp $
25  * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_poisson_pseudo_1d/ope_poisson_pseudo_1d_solp.C,v 1.4 2016/12/05 16:18:13 j_novak Exp $
26  *
27  */
28 #include <cmath>
29 #include <cstdlib>
30 
31 #include "proto.h"
32 #include "ope_elementary.h"
33 //--------------------------------------------------
34 // Version Tbl --> Tbl a 1D pour la source
35 //--------------------------------------------------
36 
37 
38 namespace Lorene {
39 Tbl _cl_poisson_pseudo_1d_pas_prevu (const Tbl &source) {
40  cout << "Combinaison lineaire pas prevue..." << endl ;
41  abort() ;
42  exit(-1) ;
43  return source;
44 }
45 
46 
47 
48  //-------------------
49  //-- R_CHEB -------
50  //--------------------
51 
52 Tbl _cl_poisson_pseudo_1d_r_cheb (const Tbl &source) {
53  Tbl barre(source) ;
54  int n = source.get_dim(0) ;
55 
56  int dirac = 1 ;
57  for (int i=0 ; i<n-2 ; i++) {
58  barre.set(i) = ((1+dirac)*source(i)-source(i+2))
59  /(i+1) ;
60  if (i==0) dirac = 0 ;
61  }
62 
63  Tbl res(barre) ;
64  for (int i=0 ; i<n-4 ; i++)
65  res.set(i) = barre(i)-barre(i+2) ;
66  return res ;
67 }
68 
69  //-------------------
70  //-- R_CHEBP -----
71  //-------------------
72 
73 Tbl _cl_poisson_pseudo_1d_r_chebp (const Tbl &source) {
74  Tbl barre(source) ;
75  int n = source.get_dim(0) ;
76 
77  int dirac = 1 ;
78  for (int i=0 ; i<n-2 ; i++) {
79  barre.set(i) = (1+dirac)*source(i)-source(i+2) ;
80  if (i==0) dirac = 0 ;
81  }
82 
83  Tbl tilde(barre) ;
84  for (int i=0 ; i<n-4 ; i++)
85  tilde.set(i) = barre(i)-barre(i+2) ;
86 
87  Tbl res(tilde) ;
88  for (int i=0 ; i<n-4 ; i++)
89  res.set(i) = tilde(i)-tilde(i+1) ;
90 
91  return res ;
92 }
93 
94 
95  //-------------------
96  //-- R_CHEBI -----
97  //-------------------
98 
99 Tbl _cl_poisson_pseudo_1d_r_chebi (const Tbl &source) {
100  Tbl barre(source) ;
101  int n = source.get_dim(0) ;
102 
103  for (int i=0 ; i<n-2 ; i++)
104  barre.set(i) = source(i)-source(i+2) ;
105 
106  Tbl tilde(barre) ;
107  for (int i=0 ; i<n-4 ; i++)
108  tilde.set(i) = barre(i)-barre(i+2) ;
109 
110  Tbl res(tilde) ;
111  for (int i=0 ; i<n-4 ; i++)
112  res.set(i) = tilde(i)-tilde(i+1) ;
113 
114  return res ;
115 }
116 
117 
118 
119 
120  //----------------------------
121  //- Routine a appeler ---
122  //------------------------------
123 
124 Tbl cl_poisson_pseudo_1d (const Tbl &source, int base_r) {
125 
126  // Routines de derivation
127  static Tbl (*cl_poisson_pseudo_1d[MAX_BASE])(const Tbl &) ;
128  static int nap = 0 ;
129 
130  // Premier appel
131  if (nap==0) {
132  nap = 1 ;
133  for (int i=0 ; i<MAX_BASE ; i++) {
134  cl_poisson_pseudo_1d[i] = _cl_poisson_pseudo_1d_pas_prevu ;
135  }
136  // Les routines existantes
137  cl_poisson_pseudo_1d[R_CHEB >> TRA_R] = _cl_poisson_pseudo_1d_r_cheb ;
138  cl_poisson_pseudo_1d[R_CHEBP >> TRA_R] = _cl_poisson_pseudo_1d_r_chebp ;
139  cl_poisson_pseudo_1d[R_CHEBI >> TRA_R] = _cl_poisson_pseudo_1d_r_chebi ;
140  }
141 
142  Tbl res(cl_poisson_pseudo_1d[base_r](source)) ;
143  return res ;
144 }
145 
146 
147  //------------------------------------
148  // Routine pour les cas non prevus --
149  //------------------------------------
150 Tbl _solp_poisson_pseudo_1d_pas_prevu (const Matrice &, const Matrice &,
151  double, double, const Tbl &) {
152  cout << " Solution homogene pas prevue ..... : "<< endl ;
153  abort() ;
154  exit(-1) ;
155  Tbl res(1) ;
156  return res;
157 }
158 
159 
160  //-------------------
161  //-- R_CHEB ------
162  //-------------------
163 
164 Tbl _solp_poisson_pseudo_1d_r_cheb (const Matrice &lap,
165  const Matrice &nondege,
166  double alpha, double beta,
167  const Tbl &source) {
168 
169  int n = lap.get_dim(0) ;
170  int dege = n-nondege.get_dim(0) ;
171  assert (dege ==2) ;
172 
173  Tbl source_aux(source*alpha*alpha) ;
174  Tbl xso(source_aux) ;
175  Tbl xxso(source_aux) ;
176  multx_1d(n, &xso.t, R_CHEB) ;
177  multx_1d(n, &xxso.t, R_CHEB) ;
178  multx_1d(n, &xxso.t, R_CHEB) ;
179  source_aux = beta*beta/alpha/alpha*source_aux+2*beta/alpha*xso+xxso ;
180  source_aux = cl_poisson_pseudo_1d(source_aux, R_CHEB) ;
181 
182  Tbl so(n-dege) ;
183  so.set_etat_qcq() ;
184  for (int i=0 ; i<n-dege ; i++)
185  so.set(i) = source_aux(i) ;
186 
187  Tbl auxi(nondege.inverse(so)) ;
188 
189  Tbl res(n) ;
190  res.set_etat_qcq() ;
191  for (int i=dege ; i<n ; i++)
192  res.set(i) = auxi(i-dege) ;
193 
194  for (int i=0 ; i<dege ; i++)
195  res.set(i) = 0 ;
196  return res ;
197 }
198 
199 
200  //-------------------
201  //-- R_CHEBP -----
202  //-------------------
203 
204 Tbl _solp_poisson_pseudo_1d_r_chebp (const Matrice &lap,
205  const Matrice &nondege,
206  double alpha, double , const Tbl &source) {
207 
208  int n = lap.get_dim(0) ;
209  int dege = n-nondege.get_dim(0) ;
210  assert ((dege==2) || (dege == 1)) ;
211  Tbl source_aux(alpha*alpha*source) ;
212  source_aux = cl_poisson_pseudo_1d(source_aux, R_CHEBP) ;
213 
214  Tbl so(n-dege) ;
215  so.set_etat_qcq() ;
216  for (int i=0 ; i<n-dege ; i++)
217  so.set(i) = source_aux(i);
218 
219  Tbl auxi(nondege.inverse(so)) ;
220 
221  Tbl res(n) ;
222  res.set_etat_qcq() ;
223  for (int i=dege ; i<n ; i++)
224  res.set(i) = auxi(i-dege) ;
225 
226  for (int i=0 ; i<dege ; i++)
227  res.set(i) = 0 ;
228 
229  if (dege==2) {
230  double somme = 0 ;
231  for (int i=0 ; i<n ; i++)
232  if (i%2 == 0)
233  somme -= res(i) ;
234  else somme += res(i) ;
235  res.set(0) = somme ;
236  return res ;
237  }
238  else return res ;
239 }
240 
241 
242  //-------------------
243  //-- R_CHEBI -----
244  //-------------------
245 
246 Tbl _solp_poisson_pseudo_1d_r_chebi (const Matrice &lap,
247  const Matrice &nondege,
248  double alpha, double, const Tbl &source) {
249 
250 
251  int n = lap.get_dim(0) ;
252  int dege = n-nondege.get_dim(0) ;
253  assert ((dege==2) || (dege == 1)) ;
254  Tbl source_aux(source*alpha*alpha) ;
255  source_aux = cl_poisson_pseudo_1d(source_aux, R_CHEBI) ;
256 
257  Tbl so(n-dege) ;
258  so.set_etat_qcq() ;
259  for (int i=0 ; i<n-dege ; i++)
260  so.set(i) = source_aux(i);
261 
262  Tbl auxi(nondege.inverse(so)) ;
263 
264  Tbl res(n) ;
265  res.set_etat_qcq() ;
266  for (int i=dege ; i<n ; i++)
267  res.set(i) = auxi(i-dege) ;
268 
269  for (int i=0 ; i<dege ; i++)
270  res.set(i) = 0 ;
271 
272  if (dege==2) {
273  double somme = 0 ;
274  for (int i=0 ; i<n ; i++)
275  if (i%2 == 0)
276  somme -= (2*i+1)*res(i) ;
277  else somme += (2*i+1)*res(i) ;
278  res.set(0) = somme ;
279  return res ;
280  }
281  else return res ;
282 }
283 
284 
286 
287  if (non_dege == 0x0)
288  do_non_dege() ;
289 
290  // Routines de derivation
291  static Tbl (*solp_poisson_pseudo_1d[MAX_BASE]) (const Matrice&,
292  const Matrice&,
293  double, double,const Tbl&) ;
294  static int nap = 0 ;
295 
296  // Premier appel
297  if (nap==0) {
298  nap = 1 ;
299  for (int i=0 ; i<MAX_BASE ; i++) {
300  solp_poisson_pseudo_1d[i] = _solp_poisson_pseudo_1d_pas_prevu ;
301  }
302  // Les routines existantes
303  solp_poisson_pseudo_1d[R_CHEB >> TRA_R] = _solp_poisson_pseudo_1d_r_cheb ;
304  solp_poisson_pseudo_1d[R_CHEBP >> TRA_R] = _solp_poisson_pseudo_1d_r_chebp ;
305  solp_poisson_pseudo_1d[R_CHEBI >> TRA_R] = _solp_poisson_pseudo_1d_r_chebi ;
306  }
307 
308  Tbl res(solp_poisson_pseudo_1d[base_r] (*ope_mat, *non_dege,
309  alpha, beta, so)) ;
310 
311  Tbl valeurs (val_solp (res, alpha, base_r)) ;
312  valeurs *= sqrt(double(2)) ;
313  sp_plus = valeurs(0) ;
314  sp_minus = valeurs(1) ;
315  dsp_plus = valeurs(2) ;
316  dsp_minus = valeurs(3) ;
317 
318  return res ;
319 }
320 }
double alpha
Parameter of the associated mapping.
virtual void do_non_dege() const
Computes the non-degenerated matrix of the operator.
virtual Tbl get_solp(const Tbl &so) const
Computes the particular solution, given the source so .
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
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.
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.
#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.
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
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
Basic array class.
Definition: tbl.h:164
#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