LORENE
ope_helmholtz_minus_2d_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_helmholtz_minus_2d_solp.C,v 1.4 2016/12/05 16:18:11 j_novak Exp $
25  * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_helmholtz_minus_2d/ope_helmholtz_minus_2d_solp.C,v 1.4 2016/12/05 16:18:11 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_helmholtz_minus_2d_pas_prevu (const Tbl & source, int) {
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_helmholtz_minus_2d_r_cheb (const Tbl &source, int) {
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 
71  //-------------------
72  //-- R_CHEBU -----
73  //-------------------
74 Tbl _cl_helmholtz_minus_2d_r_chebu_deux(const Tbl&) ;
75 
76 Tbl _cl_helmholtz_minus_2d_r_chebu (const Tbl &source, int puis) {
77 
78  int n=source.get_dim(0) ;
79  Tbl res(n) ;
80  res.set_etat_qcq() ;
81 
82  switch(puis) {
83  case 2 :
84  res = _cl_helmholtz_minus_2d_r_chebu_deux(source) ;
85  break ;
86 
87  default :
88  abort() ;
89  exit(-1) ;
90  }
91  return res ;
92 }
93 
94 // Cas dzpuis = 2 ;
95 Tbl _cl_helmholtz_minus_2d_r_chebu_deux (const Tbl &source) {
96 
97  Tbl barre(source) ;
98  int n = source.get_dim(0) ;
99 
100  int dirac = 1 ;
101  for (int i=0 ; i<n-2 ; i++) {
102  barre.set(i) = ((1+dirac)*source(i)-source(i+2)) ;
103  if (i==0) dirac = 0 ;
104  }
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 bis(tilde) ;
111  for (int i=0 ; i<n-4 ; i++)
112  bis.set(i) = (tilde(i)+tilde(i+1)) ;
113 
114  Tbl res(bis) ;
115  for (int i=0 ; i<n-4 ; i++)
116  res.set(i) = (bis(i)-bis(i+1)) ;
117 
118  return res ;
119 }
120 
121 
122  //----------------------------
123  //- Routine a appeler ---
124  //------------------------------
125 
126 Tbl cl_helmholtz_minus_2d (const Tbl &source, int puis, int base_r) {
127  // Routines de derivation
128  static Tbl (*cl_helmholtz_minus_2d[MAX_BASE])(const Tbl &, int) ;
129  static int nap = 0 ;
130 
131  // Premier appel
132  if (nap==0) {
133  nap = 1 ;
134  for (int i=0 ; i<MAX_BASE ; i++) {
135  cl_helmholtz_minus_2d[i] = _cl_helmholtz_minus_2d_pas_prevu ;
136  }
137  // Les routines existantes
138  cl_helmholtz_minus_2d[R_CHEB >> TRA_R] = _cl_helmholtz_minus_2d_r_cheb ;
139  cl_helmholtz_minus_2d[R_CHEBU >> TRA_R] = _cl_helmholtz_minus_2d_r_chebu ;
140 
141  }
142 
143  Tbl res(cl_helmholtz_minus_2d[base_r](source, puis)) ;
144  return res ;
145 }
146 
147 
148  //------------------------------------
149  // Routine pour les cas non prevus --
150  //------------------------------------
151 Tbl _solp_helmholtz_minus_2d_pas_prevu (const Matrice &, const Matrice &,
152  double, double, const Tbl &, int) {
153  cout << " Solution homogene pas prevue ..... : "<< endl ;
154  abort() ;
155  exit(-1) ;
156  Tbl res(1) ;
157  return res;
158 }
159 
160 
161  //-------------------
162  //-- R_CHEB ------
163  //-------------------
164 
165 Tbl _solp_helmholtz_minus_2d_r_cheb (const Matrice &lap, const Matrice &nondege,
166  double alpha, double beta,
167  const Tbl &source, int) {
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_helmholtz_minus_2d(source_aux, 0, 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 
197  return res ;
198 }
199 
200 
201 
202  //-------------------
203  //-- R_CHEBU -----
204  //-------------------
205 Tbl _solp_helmholtz_minus_2d_r_chebu_deux (const Matrice&, const Matrice&,
206  const Tbl&) ;
207 
208 Tbl _solp_helmholtz_minus_2d_r_chebu (const Matrice &lap, const Matrice &nondege,
209  double, double,
210  const Tbl &source, int puis) {
211  int n = lap.get_dim(0) ;
212  Tbl res(n+2) ;
213  res.set_etat_qcq() ;
214 
215  switch (puis) {
216  case 2 :
217  res = _solp_helmholtz_minus_2d_r_chebu_deux
218  (lap, nondege, source) ;
219  break ;
220  default :
221  abort() ;
222  exit(-1) ;
223  }
224 return res ;
225 }
226 
227 // Cas dzpuis = 2 ;
228 Tbl _solp_helmholtz_minus_2d_r_chebu_deux (const Matrice &lap, const Matrice &nondege,
229  const Tbl &source) {
230 
231  int n = lap.get_dim(0)+2 ;
232  int dege = n-nondege.get_dim(0) ;
233  assert (dege == 3) ;
234 
235  Tbl source_cl (cl_helmholtz_minus_2d(source, 2, R_CHEBU)) ;
236 
237  Tbl so(n-dege) ;
238  so.set_etat_qcq() ;
239  for (int i=0 ; i<n-dege ; i++)
240  so.set(i) = source_cl(i);
241 
242  Tbl sol (nondege.inverse(so)) ;
243 
244  Tbl res(n) ;
245  res.annule_hard() ;
246  for (int i=1 ; i<n-2 ; i++) {
247  res.set(i) += sol(i-1)*(2*i+3) ;
248  res.set(i+1) += -sol(i-1)*(4*i+4) ;
249  res.set(i+2) += sol(i-1)*(2*i+1) ;
250  }
251 
252  return res ;
253 }
254 
255 
257 
258  if (non_dege == 0x0)
259  do_non_dege() ;
260 
261  // Routines de derivation
262  static Tbl (*solp_helmholtz_minus_2d[MAX_BASE]) (const Matrice&, const Matrice&,
263  double, double,const Tbl&, int) ;
264  static int nap = 0 ;
265 
266  // Premier appel
267  if (nap==0) {
268  nap = 1 ;
269  for (int i=0 ; i<MAX_BASE ; i++) {
270  solp_helmholtz_minus_2d[i] = _solp_helmholtz_minus_2d_pas_prevu ;
271  }
272  // Les routines existantes
273  solp_helmholtz_minus_2d[R_CHEB >> TRA_R] = _solp_helmholtz_minus_2d_r_cheb ;
274  solp_helmholtz_minus_2d[R_CHEBU >> TRA_R] = _solp_helmholtz_minus_2d_r_chebu ;
275  }
276 
277  Tbl res(solp_helmholtz_minus_2d[base_r] (*ope_cl, *non_dege,
278  alpha, beta, so, dzpuis)) ;
279 
280  Tbl valeurs (val_solp (res, alpha, base_r)) ;
281  valeurs *= sqrt(double(2)) ;
282 
283  sp_plus = valeurs(0) ;
284  sp_minus = valeurs(1) ;
285  dsp_plus = valeurs(2) ;
286  dsp_minus = valeurs(3) ;
287 
288 
289  return res ;
290 }
291 }
double alpha
Parameter of the associated mapping.
Matrice * ope_cl
Pointer on the banded-matrix of the operator.
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
double dsp_minus
Value of the derivative of the particular solution at the inner boundary.
virtual Tbl get_solp(const Tbl &so) const
Computes the particular solution, given the source so .
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
int get_dim(int i) const
Returns the dimension of the matrix.
Definition: matrice.C:263
int dzpuis
the associated dzpuis, if in the compactified domain.
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