LORENE
ope_poisson_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_poisson_2d_solp.C,v 1.4 2016/12/05 16:18:12 j_novak Exp $
25  * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_poisson_2d/ope_poisson_2d_solp.C,v 1.4 2016/12/05 16:18:12 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_2d_pas_prevu (const Tbl &source, int puis) {
40  cout << "Combinaison lineaire pas prevue..." << endl ;
41  cout << "source : " << &source << endl ;
42  cout << "dzpuis : " << puis << endl ;
43  abort() ;
44  exit(-1) ;
45  return source;
46 }
47 
48 
49 
50  //-------------------
51  //-- R_CHEB -------
52  //--------------------
53 
54 Tbl _cl_poisson_2d_r_cheb (const Tbl &source, int) {
55  Tbl barre(source) ;
56  int n = source.get_dim(0) ;
57 
58  int dirac = 1 ;
59  for (int i=0 ; i<n-2 ; i++) {
60  barre.set(i) = ((1+dirac)*source(i)-source(i+2))
61  /(i+1) ;
62  if (i==0) dirac = 0 ;
63  }
64 
65  Tbl res(barre) ;
66  for (int i=0 ; i<n-4 ; i++)
67  res.set(i) = barre(i)-barre(i+2) ;
68  return res ;
69 }
70 
71  //-------------------
72  //-- R_CHEBP -----
73  //-------------------
74 
75 Tbl _cl_poisson_2d_r_chebp (const Tbl &source, int) {
76  Tbl barre(source) ;
77  int n = source.get_dim(0) ;
78 
79  int dirac = 1 ;
80  for (int i=0 ; i<n-2 ; i++) {
81  barre.set(i) = (1+dirac)*source(i)-source(i+2) ;
82  if (i==0) dirac = 0 ;
83  }
84 
85  Tbl tilde(barre) ;
86  for (int i=0 ; i<n-4 ; i++)
87  tilde.set(i) = barre(i)-barre(i+2) ;
88 
89  Tbl res(tilde) ;
90  for (int i=0 ; i<n-4 ; i++)
91  res.set(i) = tilde(i)-tilde(i+1) ;
92 
93  return res ;
94 }
95 
96 
97  //-------------------
98  //-- R_CHEBI -----
99  //-------------------
100 
101 Tbl _cl_poisson_2d_r_chebi (const Tbl &source, int) {
102  Tbl barre(source) ;
103  int n = source.get_dim(0) ;
104 
105  for (int i=0 ; i<n-2 ; i++)
106  barre.set(i) = source(i)-source(i+2) ;
107 
108  Tbl tilde(barre) ;
109  for (int i=0 ; i<n-4 ; i++)
110  tilde.set(i) = barre(i)-barre(i+2) ;
111 
112  Tbl res(tilde) ;
113  for (int i=0 ; i<n-4 ; i++)
114  res.set(i) = tilde(i)-tilde(i+1) ;
115 
116  return res ;
117 }
118 
119 
120  //-------------------
121  //-- R_CHEBU -----
122  //-------------------
123 Tbl _cl_poisson_2d_r_chebu_quatre(const Tbl&) ;
124 Tbl _cl_poisson_2d_r_chebu_trois(const Tbl&) ;
125 Tbl _cl_poisson_2d_r_chebu_deux(const Tbl&) ;
126 
127 Tbl _cl_poisson_2d_r_chebu (const Tbl &source, int puis) {
128 
129  int n=source.get_dim(0) ;
130  Tbl res(n) ;
131  res.set_etat_qcq() ;
132 
133  switch(puis) {
134  case 4 :
135  res = _cl_poisson_2d_r_chebu_quatre(source) ;
136  break ;
137  case 3 :
138  res = _cl_poisson_2d_r_chebu_trois (source) ;
139  break ;
140  case 2 :
141  res = _cl_poisson_2d_r_chebu_deux(source) ;
142  break ;
143 
144  default :
145  abort() ;
146  exit(-1) ;
147  }
148  return res ;
149 }
150 
151 // Cas dzpuis = 4 ;
152 Tbl _cl_poisson_2d_r_chebu_quatre (const Tbl &source) {
153  Tbl barre(source) ;
154  int n = source.get_dim(0) ;
155 
156  int dirac = 1 ;
157  for (int i=0 ; i<n-2 ; i++) {
158  barre.set(i) = ((1+dirac)*source(i)-source(i+2)) ;
159  if (i==0) dirac = 0 ;
160  }
161 
162  Tbl tilde(barre) ;
163  for (int i=0 ; i<n-4 ; i++)
164  tilde.set(i) = (barre(i)-barre(i+2)) ;
165 
166  Tbl prime(tilde) ;
167  for (int i=0 ; i<n-4 ; i++)
168  prime.set(i) = (tilde(i)-tilde(i+1)) ;
169 
170  Tbl res(prime) ;
171  for (int i=0 ; i<n-4 ; i++)
172  res.set(i) = (prime(i)-prime(i+2)) ;
173 
174  return res ;
175 }
176 // cas dzpuis = 3
177 Tbl _cl_poisson_2d_r_chebu_trois (const Tbl &source) {
178  Tbl barre(source) ;
179  int n = source.get_dim(0) ;
180 
181  int dirac = 1 ;
182  for (int i=0 ; i<n-2 ; i++) {
183  barre.set(i) = ((1+dirac)*source(i)-source(i+2)) ;
184  if (i==0) dirac = 0 ;
185  }
186 
187  Tbl tilde(barre) ;
188  for (int i=0 ; i<n-4 ; i++)
189  tilde.set(i) = (barre(i)-barre(i+2)) ;
190 
191  Tbl res(tilde) ;
192  for (int i=0 ; i<n-4 ; i++)
193  res.set(i) = (tilde(i)+tilde(i+1)) ;
194 
195  return res ;
196 }
197 
198 // Cas dzpuis = 2 ;
199 Tbl _cl_poisson_2d_r_chebu_deux (const Tbl &source) {
200  Tbl barre(source) ;
201  int n = source.get_dim(0) ;
202 
203  int dirac = 1 ;
204  for (int i=0 ; i<n-2 ; i++) {
205  barre.set(i) = ((1+dirac)*source(i)-source(i+2)) ;
206  if (i==0) dirac = 0 ;
207  }
208 
209  Tbl tilde(barre) ;
210  for (int i=0 ; i<n-4 ; i++)
211  tilde.set(i) = (barre(i)-barre(i+2)) ;
212 
213  Tbl res(tilde) ;
214  for (int i=0 ; i<n-4 ; i++)
215  res.set(i) = (tilde(i)+tilde(i+1)) ;
216  return res ;
217 }
218 
219 
220  //----------------------------
221  //- Routine a appeler ---
222  //------------------------------
223 
224 Tbl cl_poisson_2d (const Tbl &source, int puis, int base_r) {
225 
226  // Routines de derivation
227  static Tbl (*cl_poisson_2d[MAX_BASE])(const Tbl &, int) ;
228  static int nap = 0 ;
229 
230  // Premier appel
231  if (nap==0) {
232  nap = 1 ;
233  for (int i=0 ; i<MAX_BASE ; i++) {
234  cl_poisson_2d[i] = _cl_poisson_2d_pas_prevu ;
235  }
236  // Les routines existantes
237  cl_poisson_2d[R_CHEB >> TRA_R] = _cl_poisson_2d_r_cheb ;
238  cl_poisson_2d[R_CHEBU >> TRA_R] = _cl_poisson_2d_r_chebu ;
239  cl_poisson_2d[R_CHEBP >> TRA_R] = _cl_poisson_2d_r_chebp ;
240  cl_poisson_2d[R_CHEBI >> TRA_R] = _cl_poisson_2d_r_chebi ;
241  }
242 
243  Tbl res(cl_poisson_2d[base_r](source, puis)) ;
244  return res ;
245 }
246 
247 
248  //------------------------------------
249  // Routine pour les cas non prevus --
250  //------------------------------------
251 Tbl _solp_poisson_2d_pas_prevu (const Matrice &, const Matrice &,
252  double, double, const Tbl &, int) {
253  cout << " Solution homogene pas prevue ..... : "<< endl ;
254  abort() ;
255  exit(-1) ;
256  Tbl res(1) ;
257  return res;
258 }
259 
260 
261  //-------------------
262  //-- R_CHEB ------
263  //-------------------
264 
265 Tbl _solp_poisson_2d_r_cheb (const Matrice &lap, const Matrice &nondege,
266  double alpha, double beta,
267  const Tbl &source, int) {
268 
269  int n = lap.get_dim(0) ;
270  int dege = n-nondege.get_dim(0) ;
271  assert (dege ==2) ;
272 
273  Tbl source_aux(source*alpha*alpha) ;
274  Tbl xso(source_aux) ;
275  Tbl xxso(source_aux) ;
276  multx_1d(n, &xso.t, R_CHEB) ;
277  multx_1d(n, &xxso.t, R_CHEB) ;
278  multx_1d(n, &xxso.t, R_CHEB) ;
279  source_aux = beta*beta/alpha/alpha*source_aux+2*beta/alpha*xso+xxso ;
280  source_aux = cl_poisson_2d(source_aux, 0, R_CHEB) ;
281 
282  Tbl so(n-dege) ;
283  so.set_etat_qcq() ;
284  for (int i=0 ; i<n-dege ; i++)
285  so.set(i) = source_aux(i) ;
286 
287  Tbl auxi(nondege.inverse(so)) ;
288 
289  Tbl res(n) ;
290  res.set_etat_qcq() ;
291  for (int i=dege ; i<n ; i++)
292  res.set(i) = auxi(i-dege) ;
293 
294  for (int i=0 ; i<dege ; i++)
295  res.set(i) = 0 ;
296  return res ;
297 }
298 
299 
300  //-------------------
301  //-- R_CHEBP -----
302  //-------------------
303 
304 Tbl _solp_poisson_2d_r_chebp (const Matrice &lap, const Matrice &nondege,
305  double alpha, double , const Tbl &source, int) {
306 
307  int n = lap.get_dim(0) ;
308  int dege = n-nondege.get_dim(0) ;
309  assert ((dege==2) || (dege == 1)) ;
310  Tbl source_aux(alpha*alpha*source) ;
311  source_aux = cl_poisson_2d(source_aux, 0, R_CHEBP) ;
312 
313  Tbl so(n-dege) ;
314  so.set_etat_qcq() ;
315  for (int i=0 ; i<n-dege ; i++)
316  so.set(i) = source_aux(i);
317 
318  Tbl auxi(nondege.inverse(so)) ;
319 
320  Tbl res(n) ;
321  res.set_etat_qcq() ;
322  for (int i=dege ; i<n ; i++)
323  res.set(i) = auxi(i-dege) ;
324 
325  for (int i=0 ; i<dege ; i++)
326  res.set(i) = 0 ;
327 
328  if (dege==2) {
329  double somme = 0 ;
330  for (int i=0 ; i<n ; i++)
331  if (i%2 == 0)
332  somme -= res(i) ;
333  else somme += res(i) ;
334  res.set(0) = somme ;
335  return res ;
336  }
337  else return res ;
338 }
339 
340 
341  //-------------------
342  //-- R_CHEBI -----
343  //-------------------
344 
345 Tbl _solp_poisson_2d_r_chebi (const Matrice &lap, const Matrice &nondege,
346  double alpha, double, const Tbl &source, int) {
347 
348 
349  int n = lap.get_dim(0) ;
350  int dege = n-nondege.get_dim(0) ;
351  assert ((dege == 2) || (dege ==1)) ;
352  Tbl source_aux(source*alpha*alpha) ;
353  source_aux = cl_poisson_2d(source_aux, 0, R_CHEBI) ;
354 
355  Tbl so(n-dege) ;
356  so.set_etat_qcq() ;
357  for (int i=0 ; i<n-dege ; i++)
358  so.set(i) = source_aux(i);
359 
360  Tbl auxi(nondege.inverse(so)) ;
361 
362  Tbl res(n) ;
363  res.set_etat_qcq() ;
364  for (int i=dege ; i<n ; i++)
365  res.set(i) = auxi(i-dege) ;
366 
367  for (int i=0 ; i<dege ; i++)
368  res.set(i) = 0 ;
369 
370  if (dege==2) {
371  double somme = 0 ;
372  for (int i=0 ; i<n ; i++)
373  if (i%2 == 0)
374  somme -= (2*i+1)*res(i) ;
375  else somme += (2*i+1)*res(i) ;
376  res.set(0) = somme ;
377  return res ;
378  }
379  else return res ;
380 }
381 
382 
383 
384  //-------------------
385  //-- R_CHEBU -----
386  //-------------------
387 Tbl _solp_poisson_2d_r_chebu_quatre (const Matrice&, const Matrice&,
388  double, const Tbl&) ;
389 Tbl _solp_poisson_2d_r_chebu_trois (const Matrice&, const Matrice&,
390  double, const Tbl&) ;
391 Tbl _solp_poisson_2d_r_chebu_deux (const Matrice&, const Matrice&,
392  const Tbl&) ;
393 
394 Tbl _solp_poisson_2d_r_chebu (const Matrice &lap, const Matrice &nondege,
395  double alpha, double,
396  const Tbl &source, int puis) {
397  int n = lap.get_dim(0) ;
398  Tbl res(n) ;
399  res.set_etat_qcq() ;
400 
401  switch (puis) {
402 
403  case 4 :
404  res = _solp_poisson_2d_r_chebu_quatre
405  (lap, nondege, alpha, source) ;
406  break ;
407  case 3 :
408  res = _solp_poisson_2d_r_chebu_trois
409  (lap, nondege, alpha, source) ;
410  break ;
411  case 2 :
412  res = _solp_poisson_2d_r_chebu_deux
413  (lap, nondege, source) ;
414  break ;
415  default :
416  abort() ;
417  exit(-1) ;
418  }
419 return res ;
420 }
421 
422 
423 // Cas dzpuis = 4 ;
424 Tbl _solp_poisson_2d_r_chebu_quatre (const Matrice &lap, const Matrice &nondege,
425  double alpha, const Tbl &source) {
426 
427  int n = lap.get_dim(0) ;
428  int dege = n-nondege.get_dim(0) ;
429  assert ((dege==3) || (dege ==2)) ;
430  Tbl source_aux(source*alpha*alpha) ;
431  source_aux = cl_poisson_2d(source_aux, 4, R_CHEBU) ;
432 
433  Tbl so(n-dege) ;
434  so.set_etat_qcq() ;
435  for (int i=0 ; i<n-dege ; i++)
436  so.set(i) = source_aux(i);
437 
438  Tbl auxi(nondege.inverse(so)) ;
439 
440  Tbl res(n) ;
441  res.set_etat_qcq() ;
442  for (int i=dege ; i<n ; i++)
443  res.set(i) = auxi(i-dege) ;
444 
445  for (int i=0 ; i<dege ; i++)
446  res.set(i) = 0 ;
447 
448  if (dege==3) {
449  double somme = 0 ;
450  for (int i=0 ; i<n ; i++)
451  somme += i*i*res(i) ;
452  double somme_deux = somme ;
453  for (int i=0 ; i<n ; i++)
454  somme_deux -= res(i) ;
455  res.set(1) = -somme ;
456  res.set(0) = somme_deux ;
457  return res ;
458  }
459  else {
460  double somme = 0 ;
461  for (int i=0 ; i<n ; i++)
462  somme += res(i) ;
463  res.set(0) = -somme ;
464  return res ;
465  }
466 }
467 
468 // Cas dzpuis = 3 ;
469 Tbl _solp_poisson_2d_r_chebu_trois (const Matrice &lap, const Matrice &nondege,
470  double alpha, const Tbl &source) {
471 
472  int n = lap.get_dim(0) ;
473  int dege = n-nondege.get_dim(0) ;
474  assert (dege ==2) ;
475 
476  Tbl source_aux(source*alpha) ;
477  source_aux = cl_poisson_2d(source_aux, 3, R_CHEBU) ;
478 
479  Tbl so(n-dege) ;
480  so.set_etat_qcq() ;
481  for (int i=0 ; i<n-dege ; i++)
482  so.set(i) = source_aux(i);
483 
484  Tbl auxi(nondege.inverse(so)) ;
485 
486  Tbl res(n) ;
487  res.set_etat_qcq() ;
488  for (int i=dege ; i<n ; i++)
489  res.set(i) = auxi(i-dege) ;
490 
491  for (int i=0 ; i<dege ; i++)
492  res.set(i) = 0 ;
493 
494  double somme = 0 ;
495  for (int i=0 ; i<n ; i++)
496  somme += res(i) ;
497  res.set(0) = -somme ;
498  return res ;
499 }
500 
501 
502 // Cas dzpuis = 2 ;
503 Tbl _solp_poisson_2d_r_chebu_deux (const Matrice &lap, const Matrice &nondege,
504  const Tbl &source) {
505 
506  int n = lap.get_dim(0) ;
507  int dege = n-nondege.get_dim(0) ;
508  assert ((dege==1) || (dege ==2)) ;
509  Tbl source_aux(cl_poisson_2d(source, 2, R_CHEBU)) ;
510 
511  Tbl so(n-dege) ;
512  so.set_etat_qcq() ;
513  for (int i=0 ; i<n-dege ; i++)
514  so.set(i) = source_aux(i);
515 
516  Tbl auxi(nondege.inverse(so)) ;
517 
518  Tbl res(n) ;
519  res.set_etat_qcq() ;
520  for (int i=dege ; i<n ; i++)
521  res.set(i) = auxi(i-dege) ;
522 
523  for (int i=0 ; i<dege ; i++)
524  res.set(i) = 0 ;
525 
526  if (dege == 2) {
527  double somme = 0 ;
528  for (int i=0 ; i<n ; i++)
529  somme+=res(i) ;
530 
531  res.set(0) = -somme ;
532  }
533 
534  return res ;
535 }
536 
537 
538 Tbl Ope_poisson_2d::get_solp (const Tbl& so) const {
539 
540  if (non_dege == 0x0)
541  do_non_dege() ;
542 
543  // Routines de derivation
544  static Tbl (*solp_poisson_2d[MAX_BASE]) (const Matrice&, const Matrice&,
545  double, double,const Tbl&, int) ;
546  static int nap = 0 ;
547 
548  // Premier appel
549  if (nap==0) {
550  nap = 1 ;
551  for (int i=0 ; i<MAX_BASE ; i++) {
552  solp_poisson_2d[i] = _solp_poisson_2d_pas_prevu ;
553  }
554  // Les routines existantes
555  solp_poisson_2d[R_CHEB >> TRA_R] = _solp_poisson_2d_r_cheb ;
556  solp_poisson_2d[R_CHEBP >> TRA_R] = _solp_poisson_2d_r_chebp ;
557  solp_poisson_2d[R_CHEBI >> TRA_R] = _solp_poisson_2d_r_chebi ;
558  solp_poisson_2d[R_CHEBU >> TRA_R] = _solp_poisson_2d_r_chebu ;
559  }
560 
561  Tbl res(solp_poisson_2d[base_r] (*ope_mat, *non_dege,
562  alpha, beta, so, dzpuis)) ;
563 
564  Tbl valeurs (val_solp (res, alpha, base_r)) ;
565  valeurs *= sqrt(double(2)) ;
566  sp_plus = valeurs(0) ;
567  sp_minus = valeurs(1) ;
568  dsp_plus = valeurs(2) ;
569  dsp_minus = valeurs(3) ;
570 
571  return res ;
572 }
573 }
Lorene::Ope_elementary::alpha
double alpha
Parameter of the associated mapping.
Definition: ope_elementary.h:108
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
Lorene::Ope_poisson_2d::do_non_dege
virtual void do_non_dege() const
Computes the non-degenerated matrix of the operator.
Definition: ope_poisson_2d_non_dege.C:454
Lorene::Ope_elementary::beta
double beta
Parameter of the associated mapping.
Definition: ope_elementary.h:109
Lorene::Ope_poisson_2d::dzpuis
int dzpuis
the associated dzpuis, if in the compactified domain.
Definition: ope_elementary.h:727
Lorene
Lorene prototypes.
Definition: app_hor.h:67
Lorene::Ope_elementary::ope_mat
Matrice * ope_mat
Pointer on the matrix representation of the operator.
Definition: ope_elementary.h:114
Lorene::Ope_elementary::dsp_minus
double dsp_minus
Value of the derivative of the particular solution at the inner boundary.
Definition: ope_elementary.h:173
Lorene::Ope_elementary::sp_minus
double sp_minus
Value of the particular solution at the inner boundary.
Definition: ope_elementary.h:165
Lorene::Ope_elementary::base_r
int base_r
Radial basis of decomposition.
Definition: ope_elementary.h:107
Lorene::Ope_elementary::sp_plus
double sp_plus
Value of the particular solution at the outer boundary.
Definition: ope_elementary.h:169
TRA_R
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
Lorene::Ope_elementary::dsp_plus
double dsp_plus
Value of the derivative of the particular solution at the outer boundary.
Definition: ope_elementary.h:177
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Lorene::Matrice
Matrix handling.
Definition: matrice.h:152
Lorene::Tbl::get_dim
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:423
Lorene::Matrice::get_dim
int get_dim(int i) const
Returns the dimension of the matrix.
Definition: matrice.C:263
Lorene::Ope_poisson_2d::get_solp
virtual Tbl get_solp(const Tbl &so) const
Computes the particular solution, given the source so .
Definition: ope_poisson_2d_solp.C:538
Lorene::Tbl
Basic array class.
Definition: tbl.h:164
R_CHEBU
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
MAX_BASE
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
Lorene::Ope_elementary::non_dege
Matrice * non_dege
Pointer on the non-degenerated matrix of the operator.
Definition: ope_elementary.h:122
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166