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  //------------------------------------
143 namespace Lorene {
144 Tbl _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 
164 Tbl _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 
201 Tbl _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 
233 Tbl _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 
274 Tbl _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 
317 Tbl _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  }
340 return res ;
341 }
342 
343 
344 // Cas dzpuis = 4 ;
345 Tbl _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 ;
390 Tbl _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 ;
424 Tbl _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 ;
459 Tbl _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 
499 Tbl 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
Lorene prototypes.
Definition: app_hor.h:67
R_JACO02
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
TRA_R
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
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::Tbl::get_dim
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:423
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
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166