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
21char ope_poisson_2d_solp_C[] = "$Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_poisson_2d/ope_poisson_2d_solp.C,v 1.3 2014/10/13 08:53:35 j_novak Exp $" ;
22
23/*
24 * $Id: ope_poisson_2d_solp.C,v 1.3 2014/10/13 08:53:35 j_novak Exp $
25 * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_poisson_2d/ope_poisson_2d_solp.C,v 1.3 2014/10/13 08:53:35 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
38namespace Lorene {
39Tbl _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
54Tbl _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
75Tbl _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
101Tbl _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 //-------------------
123Tbl _cl_poisson_2d_r_chebu_quatre(const Tbl&) ;
124Tbl _cl_poisson_2d_r_chebu_trois(const Tbl&) ;
125Tbl _cl_poisson_2d_r_chebu_deux(const Tbl&) ;
126
127Tbl _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 ;
152Tbl _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
177Tbl _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 ;
199Tbl _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
224Tbl 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 //------------------------------------
251Tbl _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
265Tbl _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
304Tbl _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
345Tbl _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 //-------------------
387Tbl _solp_poisson_2d_r_chebu_quatre (const Matrice&, const Matrice&,
388 double, const Tbl&) ;
389Tbl _solp_poisson_2d_r_chebu_trois (const Matrice&, const Matrice&,
390 double, const Tbl&) ;
391Tbl _solp_poisson_2d_r_chebu_deux (const Matrice&, const Matrice&,
392 const Tbl&) ;
393
394Tbl _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 }
419return res ;
420}
421
422
423// Cas dzpuis = 4 ;
424Tbl _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 ;
469Tbl _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 ;
503Tbl _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
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
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}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Matrix handling.
Definition matrice.h:152
int get_dim(int i) const
Returns the dimension of the matrix.
Definition matrice.C:260
Matrice * ope_mat
Pointer on the matrix representation of the operator.
double beta
Parameter of the associated mapping.
double dsp_plus
Value of the derivative of the particular solution at the outer boundary.
double alpha
Parameter of the associated mapping.
int base_r
Radial basis of decomposition.
double sp_minus
Value of the particular solution at the inner boundary.
Matrice * non_dege
Pointer on the non-degenerated matrix of the operator.
double sp_plus
Value of the particular solution at the outer boundary.
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.
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 .
Basic array class.
Definition tbl.h:161
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition tbl.h:403
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
#define R_CHEB
base de Chebychev ordinaire (fin)
#define R_CHEBP
base de Cheb. paire (rare) seulement
Lorene prototypes.
Definition app_hor.h:64