LORENE
solp.C
1/*
2 * Copyright (c) 1999-2001 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 as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23
24
25/*
26 * $Id: solp.C,v 1.12 2016/12/05 16:18:10 j_novak Exp $
27 * $Log: solp.C,v $
28 * Revision 1.12 2016/12/05 16:18:10 j_novak
29 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
30 *
31 * Revision 1.11 2014/10/13 08:53:31 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.10 2014/10/06 15:16:10 j_novak
35 * Modified #include directives to use c++ syntax.
36 *
37 * Revision 1.9 2008/07/11 13:20:54 j_novak
38 * Miscellaneous functions for the wave equation.
39 *
40 * Revision 1.8 2008/02/18 13:53:44 j_novak
41 * Removal of special indentation instructions.
42 *
43 * Revision 1.7 2007/12/12 12:30:49 jl_cornou
44 * *** empty log message ***
45 *
46 * Revision 1.6 2004/10/05 15:44:21 j_novak
47 * Minor speed enhancements.
48 *
49 * Revision 1.5 2004/02/20 10:55:23 j_novak
50 * The versions dzpuis 5 -> 3 has been improved and polished. Should be
51 * operational now...
52 *
53 * Revision 1.4 2004/02/09 08:55:31 j_novak
54 * Corrected error in the arguments of _solp_r_chebu_cinq
55 *
56 * Revision 1.3 2004/02/06 10:53:54 j_novak
57 * New dzpuis = 5 -> dzpuis = 3 case (not ready yet).
58 *
59 * Revision 1.2 2002/10/16 14:37:12 j_novak
60 * Reorganization of #include instructions of standard C++, in order to
61 * use experimental version 3 of gcc.
62 *
63 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
64 * LORENE
65 *
66 * Revision 2.14 2000/09/07 12:54:42 phil
67 * *** empty log message ***
68 *
69 * Revision 2.13 2000/05/22 13:41:50 phil
70 * ajout du cas dzpuis == 3 ;
71 *
72 * Revision 2.12 1999/11/30 17:45:04 phil
73 * changement prototypage
74 *
75 * Revision 2.11 1999/10/12 09:44:44 phil
76 * *** empty log message ***
77 *
78 * Revision 2.10 1999/10/12 09:37:54 phil
79 * passage en const Matreice&
80 *
81 * Revision 2.9 1999/07/08 09:54:58 phil
82 * changement appel de multx_1d
83 *
84 * Revision 2.8 1999/07/07 10:05:20 phil
85 * Passage aux operateurs 1d
86 * /
87 *
88 * Revision 2.7 1999/07/02 15:04:48 phil
89 * *** empty log message ***
90 *
91 * Revision 2.6 1999/06/23 16:21:59 phil
92 * *** empty log message ***
93 *
94 * Revision 2.5 1999/06/23 12:35:06 phil
95 * ajout de dzpuis = 2
96 *
97 * Revision 2.4 1999/04/07 15:07:18 phil
98 * *** empty log message ***
99 *
100 * Revision 2.3 1999/04/07 15:06:14 phil
101 * Changement de calcul pour (-1)^n
102 *
103 * Revision 2.2 1999/04/07 14:55:46 phil
104 * Changement prototypage
105 *
106 * Revision 2.1 1999/04/07 14:36:38 phil
107 * passage par reference
108 *
109 * Revision 2.0 1999/04/07 14:11:24 phil
110 * *** empty log message ***
111 *
112 *
113 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/solp.C,v 1.12 2016/12/05 16:18:10 j_novak Exp $
114 *
115 */
116
117//fichiers includes
118#include <cstdio>
119#include <cstdlib>
120#include <cmath>
121
122#include "matrice.h"
123#include "type_parite.h"
124#include "proto.h"
125
126/*
127 * Calcul une solution particuliere a : Lap X = Y
128 *
129 * Entree :
130 * lap : matrice de l'operateur
131 * nondege : matrice non degeneree associee
132 * alpha et beta : voire mapping
133 * source : Tbl contenant les coefficients en r de Y
134 * puis : puissance dans la ZEC
135 * Sortie :
136 * Tbl contenant les coefficients de X
137 *
138 *
139 */
140 //------------------------------------
141 // Routine pour les cas non prevus --
142 //------------------------------------
143namespace Lorene {
144Tbl _solp_pas_prevu (const Matrice &lap, const Matrice &nondege, double alpha,
145 double beta, const Tbl &source, int puis) {
146 cout << " Solution particuliere pas prevue ..... : "<< endl ;
147 cout << " Laplacien : " << endl << lap << endl ;
148 cout << " Non degenere : " << endl << nondege << endl ;
149 cout << " Source : " << &source << endl ;
150 cout << " Alpha : " << alpha << endl ;
151 cout << " Beta : " << beta << endl ;
152 cout << " Puiss : " << puis << endl ;
153 abort() ;
154 exit(-1) ;
155 Tbl res(1) ;
156 return res;
157}
158
159
160 //-------------------
161 //-- R_CHEB ------
162 //-------------------
163
164Tbl _solp_r_cheb (const Matrice &lap, const Matrice &nondege, double alpha,
165 double beta, const Tbl &source, int) {
166
167 int n = lap.get_dim(0) ;
168 int dege = n-nondege.get_dim(0) ;
169 assert (dege ==2) ;
170
171 Tbl source_aux(source*alpha*alpha) ;
172 Tbl xso(source_aux) ;
173 Tbl xxso(source_aux) ;
174 multx_1d(n, &xso.t, R_CHEB) ;
175 multx_1d(n, &xxso.t, R_CHEB) ;
176 multx_1d(n, &xxso.t, R_CHEB) ;
177 source_aux = beta*beta/alpha/alpha*source_aux+2*beta/alpha*xso+xxso ;
178 source_aux = combinaison(source_aux, 0, R_CHEB) ;
179
180 Tbl so(n-dege) ;
181 so.set_etat_qcq() ;
182 for (int i=0 ; i<n-dege ; i++)
183 so.set(i) = source_aux(i) ;
184
185 Tbl auxi(nondege.inverse(so)) ;
186
187 Tbl res(n) ;
188 res.set_etat_qcq() ;
189 for (int i=dege ; i<n ; i++)
190 res.set(i) = auxi(i-dege) ;
191
192 for (int i=0 ; i<dege ; i++)
193 res.set(i) = 0 ;
194 return res ;
195}
196
197 //-------------------
198 //-- R_JACO02 ------
199 //-------------------
200
201Tbl _solp_r_jaco02 (const Matrice &lap, const Matrice &nondege, double alpha,
202 double, const Tbl &source, int) {
203
204 int n = lap.get_dim(0) ;
205 int dege = n-nondege.get_dim(0) ;
206 assert (dege ==2) ;
207
208 Tbl source_aux(source*alpha*alpha) ;
209 source_aux = combinaison(source_aux, 0, R_JACO02) ;
210
211 Tbl so(n-dege) ;
212 so.set_etat_qcq() ;
213 for (int i=0 ; i<n-dege ; i++)
214 so.set(i) = source_aux(i) ;
215
216 Tbl auxi(nondege.inverse(so)) ;
217
218 Tbl res(n) ;
219 res.set_etat_qcq() ;
220 for (int i=dege ; i<n ; i++)
221 res.set(i) = auxi(i-dege) ;
222
223 for (int i=0 ; i<dege ; i++)
224 res.set(i) = 0 ;
225 return res ;
226}
227
228
229 //-------------------
230 //-- R_CHEBP -----
231 //-------------------
232
233Tbl _solp_r_chebp (const Matrice &lap, const Matrice &nondege, double alpha,
234 double , const Tbl &source, int) {
235
236 int n = lap.get_dim(0) ;
237 int dege = n-nondege.get_dim(0) ;
238 assert ((dege==2) || (dege == 1)) ;
239 Tbl source_aux(alpha*alpha*source) ;
240 source_aux = combinaison(source_aux, 0, R_CHEBP) ;
241
242 Tbl so(n-dege) ;
243 so.set_etat_qcq() ;
244 for (int i=0 ; i<n-dege ; i++)
245 so.set(i) = source_aux(i);
246
247 Tbl auxi(nondege.inverse(so)) ;
248
249 Tbl res(n) ;
250 res.set_etat_qcq() ;
251 for (int i=dege ; i<n ; i++)
252 res.set(i) = auxi(i-dege) ;
253
254 for (int i=0 ; i<dege ; i++)
255 res.set(i) = 0 ;
256
257 if (dege==2) {
258 double somme = 0 ;
259 for (int i=0 ; i<n ; i++)
260 if (i%2 == 0)
261 somme -= res(i) ;
262 else somme += res(i) ;
263 res.set(0) = somme ;
264 return res ;
265 }
266 else return res ;
267}
268
269
270 //-------------------
271 //-- R_CHEBI -----
272 //-------------------
273
274Tbl _solp_r_chebi (const Matrice &lap, const Matrice &nondege, double alpha,
275 double, const Tbl &source, int) {
276
277
278 int n = lap.get_dim(0) ;
279 int dege = n-nondege.get_dim(0) ;
280 assert ((dege == 2) || (dege ==1)) ;
281 Tbl source_aux(source*alpha*alpha) ;
282 source_aux = combinaison(source_aux, 0, R_CHEBI) ;
283
284 Tbl so(n-dege) ;
285 so.set_etat_qcq() ;
286 for (int i=0 ; i<n-dege ; i++)
287 so.set(i) = source_aux(i);
288
289 Tbl auxi(nondege.inverse(so)) ;
290
291 Tbl res(n) ;
292 res.set_etat_qcq() ;
293 for (int i=dege ; i<n ; i++)
294 res.set(i) = auxi(i-dege) ;
295
296 for (int i=0 ; i<dege ; i++)
297 res.set(i) = 0 ;
298
299 if (dege==2) {
300 double somme = 0 ;
301 for (int i=0 ; i<n ; i++)
302 if (i%2 == 0)
303 somme -= (2*i+1)*res(i) ;
304 else somme += (2*i+1)*res(i) ;
305 res.set(0) = somme ;
306 return res ;
307 }
308 else return res ;
309}
310
311
312
313 //-------------------
314 //-- R_CHEBU -----
315 //-------------------
316
317Tbl _solp_r_chebu (const Matrice &lap, const Matrice &nondege, double alpha,
318 double, const Tbl &source, int puis) {
319 int n = lap.get_dim(0) ;
320 Tbl res(n) ;
321 res.set_etat_qcq() ;
322
323 switch (puis) {
324 case 5 :
325 res = _solp_r_chebu_cinq (lap, nondege, source) ;
326 break ;
327 case 4 :
328 res = _solp_r_chebu_quatre (lap, nondege, alpha, source) ;
329 break ;
330 case 3 :
331 res = _solp_r_chebu_trois (lap, nondege, alpha, source) ;
332 break ;
333 case 2 :
334 res = _solp_r_chebu_deux (lap, nondege, source) ;
335 break ;
336 default :
337 abort() ;
338 exit(-1) ;
339 }
340return res ;
341}
342
343
344// Cas dzpuis = 4 ;
345Tbl _solp_r_chebu_quatre (const Matrice &lap, const Matrice &nondege, double alpha,
346 const Tbl &source) {
347
348 int n = lap.get_dim(0) ;
349 int dege = n-nondege.get_dim(0) ;
350 assert ((dege==3) || (dege ==2)) ;
351 Tbl source_aux(source*alpha*alpha) ;
352 source_aux = combinaison(source_aux, 4, R_CHEBU) ;
353
354 Tbl so(n-dege) ;
355 so.set_etat_qcq() ;
356 for (int i=0 ; i<n-dege ; i++)
357 so.set(i) = source_aux(i);
358
359 Tbl auxi(nondege.inverse(so)) ;
360
361 Tbl res(n) ;
362 res.set_etat_qcq() ;
363 for (int i=dege ; i<n ; i++)
364 res.set(i) = auxi(i-dege) ;
365
366 for (int i=0 ; i<dege ; i++)
367 res.set(i) = 0 ;
368
369 if (dege==3) {
370 double somme = 0 ;
371 for (int i=0 ; i<n ; i++)
372 somme += i*i*res(i) ;
373 double somme_deux = somme ;
374 for (int i=0 ; i<n ; i++)
375 somme_deux -= res(i) ;
376 res.set(1) = -somme ;
377 res.set(0) = somme_deux ;
378 return res ;
379 }
380 else {
381 double somme = 0 ;
382 for (int i=0 ; i<n ; i++)
383 somme += res(i) ;
384 res.set(0) = -somme ;
385 return res ;
386 }
387}
388
389// Cas dzpuis = 3 ;
390Tbl _solp_r_chebu_trois (const Matrice &lap, const Matrice &nondege, double alpha,
391 const Tbl &source) {
392
393 int n = lap.get_dim(0) ;
394 int dege = n-nondege.get_dim(0) ;
395 assert (dege ==2) ;
396
397 Tbl source_aux(source*alpha) ;
398 source_aux = combinaison(source_aux, 3, R_CHEBU) ;
399
400 Tbl so(n-dege) ;
401 so.set_etat_qcq() ;
402 for (int i=0 ; i<n-dege ; i++)
403 so.set(i) = source_aux(i);
404
405 Tbl auxi(nondege.inverse(so)) ;
406
407 Tbl res(n) ;
408 res.set_etat_qcq() ;
409 for (int i=dege ; i<n ; i++)
410 res.set(i) = auxi(i-dege) ;
411
412 for (int i=0 ; i<dege ; i++)
413 res.set(i) = 0 ;
414
415 double somme = 0 ;
416 for (int i=0 ; i<n ; i++)
417 somme += res(i) ;
418 res.set(0) = -somme ;
419 return res ;
420}
421
422
423// Cas dzpuis = 2 ;
424Tbl _solp_r_chebu_deux (const Matrice &lap, const Matrice &nondege,
425 const Tbl &source) {
426
427 int n = lap.get_dim(0) ;
428 int dege = n-nondege.get_dim(0) ;
429 assert ((dege==1) || (dege ==2)) ;
430 Tbl source_aux(combinaison(source, 2, R_CHEBU)) ;
431
432 Tbl so(n-dege) ;
433 so.set_etat_qcq() ;
434 for (int i=0 ; i<n-dege ; i++)
435 so.set(i) = source_aux(i);
436
437 Tbl auxi(nondege.inverse(so)) ;
438
439 Tbl res(n) ;
440 res.set_etat_qcq() ;
441 for (int i=dege ; i<n ; i++)
442 res.set(i) = auxi(i-dege) ;
443
444 for (int i=0 ; i<dege ; i++)
445 res.set(i) = 0 ;
446
447 if (dege == 2) {
448 double somme = 0 ;
449 for (int i=0 ; i<n ; i++)
450 somme+=res(i) ;
451
452 res.set(0) = -somme ;
453 }
454
455 return res ;
456}
457
458// Cas dzpuis = 5 ;
459Tbl _solp_r_chebu_cinq (const Matrice &lap, const Matrice &nondege,
460 const Tbl &source) {
461
462 int n = lap.get_dim(0) ;
463 int dege = n-nondege.get_dim(0) ;
464
465 Tbl source_aux(combinaison(source, 5, R_CHEBU)) ;
466
467 Tbl so(n-dege) ;
468 so.set_etat_qcq() ;
469 for (int i=0 ; i<n-dege ; i++)
470 so.set(i) = source_aux(i);
471
472 Tbl auxi(nondege.inverse(so)) ;
473
474 Tbl res(n) ;
475 res.set_etat_qcq() ;
476 for (int i=dege ; i<n ; i++)
477 res.set(i) = auxi(i-dege) ;
478
479 for (int i=0 ; i<dege ; i++)
480 res.set(i) = 0 ;
481
482 if (dege == 2) {
483 double somme = 0 ;
484 for (int i=0 ; i<n ; i++)
485 somme+=res(i) ;
486
487 res.set(0) = -somme ;
488 }
489
490 return res ;
491}
492
493
494 //-------------------
495 //-- Fonction ----
496 //-------------------
497
498
499Tbl solp(const Matrice &lap, const Matrice &nondege, double alpha,
500 double beta, const Tbl &source, int puis, int base_r) {
501
502 // Routines de derivation
503 static Tbl (*solp[MAX_BASE])(const Matrice&, const Matrice&, double, double,
504 const Tbl&, int) ;
505 static int nap = 0 ;
506
507 // Premier appel
508 if (nap==0) {
509 nap = 1 ;
510 for (int i=0 ; i<MAX_BASE ; i++) {
511 solp[i] = _solp_pas_prevu ;
512 }
513 // Les routines existantes
514 solp[R_CHEB >> TRA_R] = _solp_r_cheb ;
515 solp[R_CHEBU >> TRA_R] = _solp_r_chebu ;
516 solp[R_CHEBP >> TRA_R] = _solp_r_chebp ;
517 solp[R_CHEBI >> TRA_R] = _solp_r_chebi ;
518 solp[R_JACO02 >> TRA_R] = _solp_r_jaco02 ;
519 }
520
521 return solp[base_r](lap, nondege, alpha, beta, source, puis) ;
522}
523}
Lorene::Matrice
Matrix handling.
Definition matrice.h:152
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::get_dim
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i]).
Definition tbl.h:403
MAX_BASE
#define MAX_BASE
Nombre max. de bases differentes.
Definition type_parite.h:144
R_CHEBU
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition type_parite.h:180
R_JACO02
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition type_parite.h:188
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition type_parite.h:170
TRA_R
#define TRA_R
Translation en R, used for a bitwise shift (in hex).
Definition type_parite.h:158
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition type_parite.h:166
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
Lorene
Lorene prototypes.
Definition app_hor.h:67