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 char 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  //------------------------------------
140 namespace Lorene {
141 Tbl _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 
161 Tbl _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 
198 Tbl _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 
230 Tbl _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 
271 Tbl _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 
314 Tbl _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  }
337 return res ;
338 }
339 
340 
341 // Cas dzpuis = 4 ;
342 Tbl _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 ;
387 Tbl _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 ;
421 Tbl _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 ;
456 Tbl _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 
496 Tbl 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 }
Lorene prototypes.
Definition: app_hor.h:64
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:403
#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
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166