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