LORENE
poisson_compact.C
1 /*
2  * Copyright (c) 2000-2006 Philippe Grandclement
3  * Copyright (c) 2007 Michal Bejger
4  * Copyright (c) 2007 Eric Gourgoulhon
5  *
6  * This file is part of LORENE.
7  *
8  * LORENE is free software; you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation; either version 2 of the License, or
11  * (at your option) any later version.
12  *
13  * LORENE is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with LORENE; if not, write to the Free Software
20  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21  *
22  */
23 
24 
25 char poisson_compact_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_compact.C,v 1.6 2014/10/13 08:53:29 j_novak Exp $" ;
26 
27 /*
28  * $Id: poisson_compact.C,v 1.6 2014/10/13 08:53:29 j_novak Exp $
29  * $Log: poisson_compact.C,v $
30  * Revision 1.6 2014/10/13 08:53:29 j_novak
31  * Lorene classes and functions now belong to the namespace Lorene.
32  *
33  * Revision 1.5 2014/10/06 15:16:09 j_novak
34  * Modified #include directives to use c++ syntax.
35  *
36  * Revision 1.4 2007/10/16 21:54:23 e_gourgoulhon
37  * Added new function sol_poisson_compact (multi-domain version).
38  *
39  * Revision 1.3 2006/09/05 13:39:46 p_grandclement
40  * update of the bin_ns_bh project
41  *
42  * Revision 1.2 2002/10/16 14:37:12 j_novak
43  * Reorganization of #include instructions of standard C++, in order to
44  * use experimental version 3 of gcc.
45  *
46  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
47  * LORENE
48  *
49  * Revision 2.2 2000/03/16 16:28:06 phil
50  * Version entirement revue et corrigee
51  *
52  * Revision 2.1 2000/03/09 13:51:55 phil
53  * *** empty log message ***
54  *
55  * Revision 2.0 2000/03/09 13:44:56 phil
56  * *** empty log message ***
57  *
58  *
59  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_compact.C,v 1.6 2014/10/13 08:53:29 j_novak Exp $
60  *
61  */
62 
63 // Headers C
64 #include <cstdlib>
65 #include <cmath>
66 #include <cassert>
67 
68 // Headers Lorene
69 #include "map.h"
70 #include "diff.h"
71 #include "matrice.h"
72 #include "type_parite.h"
73 #include "proto.h"
74 #include "base_val.h"
75 #include "utilitaires.h"
76 
78  // Single domain version //
80 
81 /*
82  * Cette fonction resout, dans le noyau :
83  * a*(1-xi^2)*lap(uu)+b*xi*duu/dxi = source
84  * avec a>0 et b<0 ;
85  * Pour le stokage des operateurs, il faut faire reamorce = true au
86  * debut d'un nouveau calcul.
87  */
88 
89 namespace Lorene {
90 Mtbl_cf sol_poisson_compact(const Mtbl_cf& source, double a, double b,
91  bool reamorce) {
92 
93  // Verifications :
94  assert (source.get_etat() != ETATNONDEF) ;
95 
96  assert (a>0) ;
97  assert (b<0) ;
98 
99  // Les tableaux de stockage :
100  const int nmax = 200 ;
101  static Matrice* tab_op[nmax] ;
102  static int nb_deja_fait = 0 ;
103  static int l_deja_fait[nmax] ;
104  static int n_deja_fait[nmax] ;
105 
106  if (reamorce) {
107  for (int i=0 ; i<nb_deja_fait ; i++)
108  delete tab_op[i] ;
109  nb_deja_fait = 0 ;
110  }
111 
112  int nz = source.get_mg()->get_nzone() ;
113 
114  // Pour le confort (on ne travaille que dans le noyau) :
115  int nr = source.get_mg()->get_nr(0) ;
116  int nt = source.get_mg()->get_nt(0) ;
117  int np = source.get_mg()->get_np(0) ;
118 
119  int l_quant ;
120  int m_quant ;
121  int base_r ;
122 
123  // La solution ...
124  Mtbl_cf solution(source.get_mg(), source.base) ;
125  solution.set_etat_qcq() ;
126  solution.t[0]->set_etat_qcq() ;
127 
128  for (int k=0 ; k<np+1 ; k++)
129  for (int j=0 ; j<nt ; j++)
130  if (nullite_plm(j, nt, k, np, source.base) == 1)
131  {
132  // calcul des nombres quantiques :
133  donne_lm(nz, 0, j, k, source.base, m_quant, l_quant, base_r) ;
134 
135 
136  //On gere le cas l_quant == 0 (c'est bien simple y'en a pas !)
137  if (l_quant == 0) {
138  for (int i=0 ; i<nr ; i++)
139  solution.set(0, k, j, i) = 0 ;
140  }
141 
142  // cas l_quant != 0
143  else {
144  // On determine si la matrice a deja ete calculee :
145  int indice = -1 ;
146 
147  // Le cas l==1 est non singulier : pas de base de Gelerkin
148  int taille = (l_quant == 1) ? nr : nr-1 ;
149 
150  Matrice operateur (taille, taille) ;
151  for (int conte=0 ; conte<nb_deja_fait ; conte++)
152  if ((l_deja_fait[conte]== l_quant) && (n_deja_fait[conte] == nr))
153  indice = conte ;
154 
155  if (indice == -1) {
156  if (nb_deja_fait >= nmax) {
157  cout << "sol_poisson_compact : trop de matrices ..." << endl;
158  abort() ;
159  }
160  // Calcul a faire :
161  operateur = a*lap_cpt_mat (nr, l_quant, base_r)
162  + b*xdsdx_mat(nr, l_quant, base_r) ;
163  operateur = combinaison_cpt (operateur, l_quant, base_r) ;
164 
165  l_deja_fait[nb_deja_fait] = l_quant ;
166  n_deja_fait[nb_deja_fait] = nr ;
167  tab_op[nb_deja_fait] = new Matrice(operateur) ;
168 
169  nb_deja_fait++ ;
170  }
171  else {
172  // rien a faire :
173  operateur = *tab_op[indice] ;
174  }
175 
176  // La source :
177  Tbl so(taille) ;
178  so.set_etat_qcq() ;
179  for (int i=0 ; i<taille ; i++)
180  so.set(i) = source(0, k, j, i) ;
181  so = combinaison_cpt (so, base_r) ;
182 
183  Tbl part (operateur.inverse(so)) ;
184 
185  if (taille == nr)
186  for (int i=0 ; i<nr ; i++)
187  solution.set(0, k, j, i) = part(i) ; // cas l==1
188  else {
189  solution.set(0, k, j, nr-1) = 0 ;
190  for (int i=nr-2 ; i>=0 ; i--)
191  if (base_r == R_CHEBP) { //Gelerkin pair
192  solution.set(0, k, j, i) = part(i) ;
193  solution.set(0, k, j, i+1) += part(i) ;
194  }
195  else { //Gelerkin impaire
196  solution.set(0, k, j, i) = part(i)*(2*i+3) ;
197  solution.set(0, k, j, i+1) += part(i)*(2*i+1) ;
198  }
199  }
200  }
201  }
202  else // cas ou nullite_plm = 0 :
203  for (int i=0 ; i<nr ; i++)
204  solution.set(0, k, j, i) = 0 ;
205 
206  // Mise a zero du coefficient (inusite) k=np+1
207  for (int j=0; j<nt; j++)
208  for(int i=0 ; i<nr ; i++)
209  solution.set(0, np+1, j, i) = 0 ;
210 
211  for (int zone=1 ; zone<nz ; zone++)
212  solution.t[zone]->set_etat_zero() ;
213 
214  return solution ;
215 }
216 
217 
218 
220  // Multi domain version //
222 
223 
224 Mtbl_cf sol_poisson_compact(const Map_af& mapping, const Mtbl_cf& source,
225  const Tbl& ac, const Tbl& bc, bool ) {
226 
227  // Number of domains inside the star :
228  int nzet = ac.get_dim(1) ;
229  assert(nzet<=source.get_mg()->get_nzone()) ;
230 
231  // Some checks
232  assert (source.get_mg()->get_type_r(0) == RARE) ;
233  for (int l=1 ; l<nzet ; l++)
234  assert(source.get_mg()->get_type_r(l) == FIN) ;
235 
236  // Spectral bases
237  const Base_val& base = source.base ;
238 
239  // Result
240  Mtbl_cf resultat(source.get_mg(), base) ;
241  resultat.annule_hard() ;
242 
243  // donnees sur la zone
244  int nr, nt, np ;
245  int base_r ;
246  double alpha, beta, echelle ;
247  int l_quant, m_quant;
248 
249  // Determination of the size of the systeme :
250  int size = 0 ;
251  int max_nr = 0 ;
252  for (int l=0 ; l<nzet ; l++) {
253  nr = mapping.get_mg()->get_nr(l) ;
254  size += nr ;
255  if (nr > max_nr) max_nr = nr ;
256  }
257 
258  Matrice systeme (size, size) ;
259  systeme.set_etat_qcq() ;
260  Tbl sec_membre (size) ;
261  Tbl soluce (size) ;
262  soluce.set_etat_qcq() ;
263 
264  np = mapping.get_mg()->get_np(0) ;
265  nt = mapping.get_mg()->get_nt(0) ;
266  Matrice* work ;
267 
268  // On bosse pour chaque l, m :
269  for (int k=0 ; k<np+1 ; k++)
270  for (int j=0 ; j<nt ; j++)
271  if (nullite_plm(j, nt, k, np, base) == 1) {
272 
273  systeme.annule_hard() ;
274  sec_membre.annule_hard() ;
275 
276  int column_courant = 0 ;
277  int ligne_courant = 0 ;
278 
279  //--------------------------
280  // NUCLEUS
281  //--------------------------
282 
283  nr = mapping.get_mg()->get_nr(0) ;
284  alpha = mapping.get_alpha()[0] ;
285  base.give_quant_numbers (0, k, j, m_quant, l_quant, base_r) ;
286 
287  if (l_quant == 0) {
288  for (int i=0 ; i<size ; i++)
289  soluce.set(i) = 0 ;
290  }
291  else {
292 
293  Diff_dsdx2 d2_n(base_r, nr) ; // suffix _n stands for "nucleus"
294  Diff_sxdsdx sxd_n(base_r, nr) ;
295  Diff_sx2 sx2_n(base_r, nr) ;
296  Diff_x2dsdx2 x2d2_n(base_r,nr) ;
297  Diff_xdsdx xd_n(base_r,nr) ;
298  Diff_id id_n(base_r,nr) ;
299 
300  work = new Matrice( ac(0,0) * ( d2_n + 2.*sxd_n -l_quant*(l_quant+1)*sx2_n )
301  + ac(0,2) * ( x2d2_n + 2.*xd_n -l_quant*(l_quant+1)*id_n )
302  + alpha * bc(0,1) * xd_n ) ;
303 
304  // cout << *work << endl ;
305  // arrete() ;
306 
307  // regularity conditions :
308  int nbr_cl = 0 ;
309  if (l_quant > 1) {
310  nbr_cl = 1 ;
311  if (l_quant%2==0) {
312  //Even case
313  for (int col=0 ; col<nr ; col++)
314  if (col%2==0)
315  systeme.set(ligne_courant, col+column_courant) = 1 ;
316  else
317  systeme.set(ligne_courant, col+column_courant) = -1 ;
318  }
319  else {
320  //Odd case
321  for (int col=0 ; col<nr ; col++)
322  if (col%2==0)
323  systeme.set(ligne_courant, col+column_courant) = 2*col+1 ;
324  else
325  systeme.set(ligne_courant, col+column_courant) = -(2*col+1) ;
326  }
327  ligne_courant ++ ;
328  }
329 
330  // L'operateur :
331  for (int lig=0 ; lig<nr-1-nbr_cl ; lig++) {
332  for (int col=0 ; col<nr ; col++)
333  systeme.set(lig+ligne_courant,col+column_courant) = (*work)(lig,col) ;
334  sec_membre.set(lig+ligne_courant) = alpha*alpha*source(0, k, j, lig) ;
335  }
336 
337  delete work ;
338  ligne_courant += nr-1-nbr_cl ;
339 
340  // Le raccord :
341  for (int col=0 ; col<nr ; col++) {
342  // La fonction
343  systeme.set(ligne_courant, col+column_courant) = 1 ;
344  // Sa dérivée :
345  if (l_quant%2==0)
346  systeme.set(ligne_courant+1, col+column_courant) = 4*col*col/alpha ;
347  else
348  systeme.set(ligne_courant+1, col+column_courant) = (2*col+1)*(2*col+1)/alpha ;
349  }
350  column_courant += nr ;
351 
352  //--------------------------
353  // SHELLS
354  //--------------------------
355 
356  for (int l=1 ; l<nzet ; l++) {
357 
358  nr = mapping.get_mg()->get_nr(l) ;
359  alpha = mapping.get_alpha()[l] ;
360  beta = mapping.get_beta()[l] ;
361  echelle = beta/alpha ;
362  double bsa = echelle ;
363  double bsa2 = bsa*bsa ;
364 
365  base.give_quant_numbers (l, k, j, m_quant, l_quant, base_r) ;
366 
367  Diff_dsdx dx(base_r, nr) ;
368  Diff_xdsdx xdx(base_r, nr) ;
369  Diff_x2dsdx x2dx(base_r, nr) ;
370  Diff_x3dsdx x3dx(base_r, nr) ;
371  Diff_dsdx2 dx2(base_r, nr) ;
372  Diff_xdsdx2 xdx2(base_r, nr) ;
373  Diff_x2dsdx2 x2dx2(base_r, nr) ;
374  Diff_x3dsdx2 x3dx2(base_r, nr) ;
375  Diff_id id(base_r,nr) ;
376  Diff_mx mx(base_r,nr) ;
377 
378  work = new Matrice ( ac(l,0) * ( x2dx2 + 2.*bsa*xdx2 + bsa2*dx2 + 2.*xdx
379  + 2.*bsa*dx - l_quant*(l_quant+1)*id )
380  + ac(l,1) * ( x3dx2 + 2.*bsa*x2dx2 + bsa2*xdx2 + 2.*x2dx
381  + 2.*bsa*xdx - l_quant*(l_quant+1) *mx )
382  + alpha * ( bc(l,0) * ( x2dx + 2.*bsa*xdx + bsa2*dx )
383  + bc(l,1) * ( x3dx + 2.*bsa*x2dx + bsa2*xdx ) ) ) ;
384 
385  // matching with previous domain :
386  for (int col=0 ; col<nr ; col++) {
387  // La fonction
388  if (col%2==0)
389  systeme.set(ligne_courant, col+column_courant) = -1 ;
390  else
391  systeme.set(ligne_courant, col+column_courant) = 1 ;
392  // Sa dérivée :
393  if (col%2==0)
394  systeme.set(ligne_courant+1, col+column_courant) = col*col/alpha ;
395  else
396  systeme.set(ligne_courant+1, col+column_courant) = -col*col/alpha ;
397  }
398  ligne_courant += 2 ;
399 
400  // source must be multiplied by (x+echelle)^2
401  Tbl source_aux(nr) ;
402  source_aux.set_etat_qcq() ;
403  for (int i=0 ; i<nr ; i++)
404  source_aux.set(i) = source(l,k,j,i)*alpha*alpha ;
405  Tbl xso(source_aux) ;
406  Tbl xxso(source_aux) ;
407  multx_1d(nr, &xso.t, R_CHEB) ;
408  multx_1d(nr, &xxso.t, R_CHEB) ;
409  multx_1d(nr, &xxso.t, R_CHEB) ;
410  source_aux = beta*beta/alpha/alpha*source_aux+2*beta/alpha*xso+xxso ;
411 
412  // L'operateur :
413 
414  for (int lig=0 ; lig<nr-1 ; lig++) {
415  for (int col=0 ; col<nr ; col++)
416  systeme.set(lig+ligne_courant,col+column_courant) = (*work)(lig,col) ;
417  sec_membre.set(lig+ligne_courant) = source_aux(lig) ;
418  }
419 
420  // cout << *work << endl ;
421  // arrete() ;
422 
423  delete work ;
424  ligne_courant += nr-2 ;
425 
426  if (l<nzet-1) { // if this not the last shell
427  // matching with the next domain
428  for (int col=0 ; col<nr ; col++) {
429  // La fonction
430  systeme.set(ligne_courant, col+column_courant) = 1 ;
431  // Sa dérivée :
432  systeme.set(ligne_courant+1, col+column_courant) = col*col/alpha ;
433  }
434  }
435 
436  column_courant += nr ;
437 
438  } // end loop on the shells (index l)
439 
440  // cout << "systeme : " << systeme << endl ;
441  // arrete() ;
442 
443  // Solving the system:
444 
445  systeme.set_band (size, size) ;
446  systeme.set_lu() ;
447 
448  // cout << "Determinant: " << systeme.determinant() << endl ;
449  // cout << "Eigenvalues: " << systeme.val_propre() << endl ;
450 
451  soluce = systeme.inverse(sec_membre) ;
452 
453  } // end case l_quant != 0
454 
455  // cout << "soluce: " << soluce << endl ;
456  // arrete() ;
457 
458  // On range :
459  int conte = 0 ;
460  for (int l=0 ; l<nzet ; l++) {
461  nr = mapping.get_mg()->get_nr(l) ;
462  for (int i=0 ; i<nr ; i++) {
463  resultat.set(l,k,j,i) = soluce(conte) ;
464  conte++ ;
465  }
466  }
467 
468  } // end case nullite_plm == 1
469 
470  return resultat ;
471 
472 }
473 
474 
475 
476 
477 
478 
479 
480 
481 
482 
483 
484 
485 
486 
487 
488 }
Lorene prototypes.
Definition: app_hor.h:64
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Base_val base
Bases of the spectral expansions.
Definition: mtbl_cf.h:200
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166