LORENE
poisson_tau.C
1 /*
2  * Copyright (c) 2005 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 poisson_tau_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_tau.C,v 1.10 2014/10/13 08:53:30 j_novak Exp $" ;
24 
25 /*
26  * $Id: poisson_tau.C,v 1.10 2014/10/13 08:53:30 j_novak Exp $
27  * $Log: poisson_tau.C,v $
28  * Revision 1.10 2014/10/13 08:53:30 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.9 2014/10/06 15:16:09 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.8 2013/06/05 15:10:43 j_novak
35  * Suppression of FINJAC sampling in r. This Jacobi(0,2) base is now
36  * available by setting colloc_r to BASE_JAC02 in the Mg3d constructor.
37  *
38  * Revision 1.7 2008/08/27 08:51:15 jl_cornou
39  * Added Jacobi(0,2) polynomials
40  *
41  * Revision 1.6 2007/12/14 10:19:34 jl_cornou
42  * *** empty log message ***
43  *
44  * Revision 1.4 2005/11/24 14:07:54 j_novak
45  * Use of Matrice::annule_hard()
46  *
47  * Revision 1.3 2005/08/26 14:02:41 p_grandclement
48  * Modification of the elliptic solver that matches with an oscillatory exterior solution
49  * small correction in Poisson tau also...
50  *
51  * Revision 1.2 2005/08/25 12:16:01 p_grandclement
52  * *** empty log message ***
53  *
54  *
55  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_tau.C,v 1.10 2014/10/13 08:53:30 j_novak Exp $
56  *
57  */
58 
59 // Header C :
60 #include <cstdlib>
61 #include <cmath>
62 
63 // Headers Lorene :
64 #include "matrice.h"
65 #include "map.h"
66 #include "proto.h"
67 #include "type_parite.h"
68 
69 
70 
71  //----------------------------------------------
72  // Version Mtbl_cf
73  //----------------------------------------------
74 
75 /*
76  *
77  * Solution de l'equation de poisson with a multi-domain Tau method
78  *
79  * Entree : mapping : le mapping affine
80  * source : les coefficients de la source qui a ete multipliee par
81  * r^4 r^3 ou r^2 dans la ZEC.
82  * La base de decomposition doit etre Ylm
83  * dzpuis : exposant de r dans le factor multiplicatif dans la ZEC
84  * Sortie : renvoie les coefficients de la solution dans la meme base de
85  * decomposition (a savoir Ylm)
86  *
87  */
88 
89 
90 namespace Lorene {
91 Mtbl_cf sol_poisson_tau(const Map_af& mapping, const Mtbl_cf& source, int dzpuis)
92 {
93 
94  // Verifications d'usage sur les zones
95  int nz = source.get_mg()->get_nzone() ;
96  assert (nz>1) ;
97  assert ((source.get_mg()->get_type_r(0) == RARE) || (source.get_mg()->get_type_r(0) == FIN)) ;
98  assert (source.get_mg()->get_type_r(nz-1) == UNSURR) ;
99  for (int l=1 ; l<nz-1 ; l++)
100  assert(source.get_mg()->get_type_r(l) == FIN) ;
101 
102  assert ((dzpuis==4) || (dzpuis==2) || (dzpuis==3)) ;
103 
104  // Bases spectrales
105  const Base_val& base = source.base ;
106 
107  // Resultat
108  Mtbl_cf resultat(source.get_mg(), base) ;
109  resultat.annule_hard() ;
110 
111  // donnees sur la zone
112  int nr, nt, np ;
113  int base_r ;
114  double alpha, beta, echelle ;
115  int l_quant, m_quant;
116 
117  // Determination of the size of the systeme :
118  int size = 0 ;
119  int max_nr = 0 ;
120  for (int l=0 ; l<nz ; l++) {
121  nr = mapping.get_mg()->get_nr(l) ;
122  size += nr ;
123  if (nr > max_nr)
124  max_nr = nr ;
125  }
126 
127  Matrice systeme (size, size) ;
128  systeme.set_etat_qcq() ;
129  Tbl sec_membre (size) ;
130 
131  np = mapping.get_mg()->get_np(0) ;
132  nt = mapping.get_mg()->get_nt(0) ;
133  Matrice* work ;
134 
135  // On bosse pour chaque l, m :
136  for (int k=0 ; k<np+1 ; k++)
137  for (int j=0 ; j<nt ; j++)
138  if (nullite_plm(j, nt, k, np, base) == 1) {
139 
140 // for (int lig=0 ; lig<size ; lig++)
141 // for (int col=0 ; col< size ; col++)
142 // systeme.set(lig,col) = 0 ;
143  systeme.annule_hard() ;
144  sec_membre.annule_hard() ;
145 
146  int column_courant = 0 ;
147  int ligne_courant = 0 ;
148 
149  //--------------------------
150  // NUCLEUS
151  //--------------------------
152 
153  nr = mapping.get_mg()->get_nr(0) ;
154  alpha = mapping.get_alpha()[0] ;
155  base.give_quant_numbers (0, k, j, m_quant, l_quant, base_r) ;
156  work = new Matrice (laplacien_mat(nr, l_quant, 0., 0, base_r)) ;
157 
158  int nbr_cl = 0 ;
159  // RARE case
160  if (source.get_mg()->get_type_r(0) == RARE) {
161  // regularity conditions :
162  if (l_quant > 1) {
163  nbr_cl = 1 ;
164  if (l_quant%2==0) {
165  //Even case
166  for (int col=0 ; col<nr ; col++)
167  if (col%2==0)
168  systeme.set(ligne_courant, col+column_courant) = 1 ;
169  else
170  systeme.set(ligne_courant, col+column_courant) = -1 ;
171  }
172  else {
173  //Odd case
174  for (int col=0 ; col<nr ; col++)
175  if (col%2==0)
176  systeme.set(ligne_courant, col+column_courant) = 2*col+1 ;
177  else
178  systeme.set(ligne_courant, col+column_courant) = -(2*col+1) ;
179  }
180  }
181  }
182 
183  // JACO02 case
184  else {
185  assert( base_r == R_JACO02) ;
186  // regularity conditions :
187  if (l_quant == 0) {
188  nbr_cl = 1 ;
189  for (int col=0 ; col<nr ; col++) {
190  systeme.set(ligne_courant, col+column_courant) = col*(col+1)*(col+2)*(col+3)/double(12)*(2*(col%2)-1);
191  }
192  }
193  else if (l_quant == 1) {
194  nbr_cl = 1 ;
195  for (int col=0 ; col<nr ; col++) {
196  systeme.set(ligne_courant, col+column_courant) = (col+1)*(col+2)/double(2)*(1-2*(col%2)) ;
197  }
198  }
199  else {
200  nbr_cl = 2 ;
201  for (int col=0 ; col<nr ; col++) {
202  systeme.set(ligne_courant, col+column_courant) = (col+1)*(col+2)/double(2)*(1-2*(col%2)) ;
203  systeme.set(ligne_courant+1, col+column_courant) = col*(col+1)*(col+2)*(col+3)/double(12)*(2*(col%2)-1) ;
204  }
205  }
206  }
207  ligne_courant += nbr_cl ;
208 
209  // L'operateur :
210  for (int lig=0 ; lig<nr-1-nbr_cl ; lig++) {
211  for (int col=0 ; col<nr ; col++)
212  systeme.set(lig+ligne_courant,col+column_courant) = (*work)(lig,col) ;
213  sec_membre.set(lig+ligne_courant) = alpha*alpha*source(0, k, j, lig) ;
214  }
215 
216  delete work ;
217  ligne_courant += nr-1-nbr_cl ;
218 
219  // Le raccord :
220  for (int col=0 ; col<nr ; col++) {
221  if (source.get_mg()->get_type_r(0) == RARE) {
222  // La fonction
223  systeme.set(ligne_courant, col+column_courant) = 1 ;
224  // Sa d�riv�e :
225  if (l_quant%2==0) {
226  systeme.set(ligne_courant+1, col+column_courant) = 4*col*col/alpha ;
227  }
228  else {
229  systeme.set(ligne_courant+1, col+column_courant) = (2*col+1)*(2*col+1)/alpha ;
230  }
231  }
232  else {
233  // La fonction
234  systeme.set(ligne_courant, col+column_courant) = 1 ;
235  // Sa dérivée :
236  systeme.set(ligne_courant+1, col+column_courant) = col*(col+3)/double(2)/alpha ;
237  }
238  }
239 
240  column_courant += nr ;
241 
242 
243 
244 
245  //--------------------------
246  // SHELLS
247  //--------------------------
248  for (int l=1 ; l<nz-1 ; l++) {
249 
250  nr = mapping.get_mg()->get_nr(l) ;
251  alpha = mapping.get_alpha()[l] ;
252  beta = mapping.get_beta()[l] ;
253  echelle = beta/alpha ;
254 
255  base.give_quant_numbers (l, k, j, m_quant, l_quant, base_r) ;
256  work = new Matrice (laplacien_mat(nr, l_quant, echelle, 0, base_r)) ;
257 
258  // matching with previous domain :
259  for (int col=0 ; col<nr ; col++) {
260  // La fonction
261  if (col%2==0)
262  systeme.set(ligne_courant, col+column_courant) = -1 ;
263  else
264  systeme.set(ligne_courant, col+column_courant) = 1 ;
265  // Sa d�riv�e :
266  if (col%2==0)
267  systeme.set(ligne_courant+1, col+column_courant) = col*col/alpha ;
268  else
269  systeme.set(ligne_courant+1, col+column_courant) = -col*col/alpha ;
270  }
271  ligne_courant += 2 ;
272 
273  // L'operateur :
274 
275  // source must be multiplied by (x+echelle)^2
276  Tbl source_aux(nr) ;
277  source_aux.set_etat_qcq() ;
278  for (int i=0 ; i<nr ; i++)
279  source_aux.set(i) = source(l,k,j,i)*alpha*alpha ;
280  Tbl xso(source_aux) ;
281  Tbl xxso(source_aux) ;
282  multx_1d(nr, &xso.t, R_CHEB) ;
283  multx_1d(nr, &xxso.t, R_CHEB) ;
284  multx_1d(nr, &xxso.t, R_CHEB) ;
285  source_aux = beta*beta/alpha/alpha*source_aux+2*beta/alpha*xso+xxso ;
286 
287  for (int lig=0 ; lig<nr-2 ; lig++) {
288  for (int col=0 ; col<nr ; col++)
289  systeme.set(lig+ligne_courant,col+column_courant) = (*work)(lig,col) ;
290  sec_membre.set(lig+ligne_courant) = source_aux(lig) ;
291  }
292 
293  delete work ;
294  ligne_courant += nr-2 ;
295  // Matching with the next domain :
296  for (int col=0 ; col<nr ; col++) {
297  // La fonction
298  systeme.set(ligne_courant, col+column_courant) = 1 ;
299  // Sa d�riv�e :
300  systeme.set(ligne_courant+1, col+column_courant) = col*col/alpha ;
301  }
302 
303  column_courant += nr ;
304  }
305 
306 
307  //--------------------------
308  // ZEC
309  //--------------------------
310  nr = mapping.get_mg()->get_nr(nz-1) ;
311  alpha = mapping.get_alpha()[nz-1] ;
312  beta = mapping.get_beta()[nz-1] ;
313 
314  base.give_quant_numbers (nz-1, k, j, m_quant, l_quant, base_r) ;
315  work = new Matrice(laplacien_mat(nr, l_quant, 0., dzpuis, base_r)) ;
316 
317  // Matching with the previous domain :
318  for (int col=0 ; col<nr ; col++) {
319  // La fonction
320  if (col%2==0)
321  systeme.set(ligne_courant, col+column_courant) = -1 ;
322  else
323  systeme.set(ligne_courant, col+column_courant) = 1 ;
324  // Sa d�riv�e :
325  if (col%2==0)
326  systeme.set(ligne_courant+1, col+column_courant) = -4*alpha*col*col ;
327  else
328  systeme.set(ligne_courant+1, col+column_courant) = 4*alpha*col*col ;
329  }
330  ligne_courant += 2 ;
331 
332  // Regularity and BC at infinity ?
333  nbr_cl =0 ;
334  switch (dzpuis) {
335  case 4 :
336  if (l_quant==0) {
337  nbr_cl = 1 ;
338  // Only BC at infinity :
339  for (int col=0 ; col<nr ; col++)
340  systeme.set(ligne_courant, col+column_courant) = 1 ;
341  }
342  else {
343  nbr_cl = 2 ;
344  // BC at infinity :
345  for (int col=0 ; col<nr ; col++)
346  systeme.set(ligne_courant, col+column_courant) = 1 ;
347  // Regularity :
348  for (int col=0 ; col<nr ; col++)
349  systeme.set(ligne_courant+1, col+column_courant) = -4*alpha*col*col ;
350  }
351  break ;
352 
353  case 3 :
354  nbr_cl = 1 ;
355  // Only BC at infinity :
356  for (int col=0 ; col<nr ; col++)
357  systeme.set(ligne_courant, col+column_courant) = 1 ;
358  break ;
359 
360  case 2 :
361  if (l_quant==0) {
362  nbr_cl = 1 ;
363  // Only BC at infinity :
364  for (int col=0 ; col<nr ; col++)
365  systeme.set(ligne_courant, col+column_courant) = 1 ;
366  }
367  break ;
368  default :
369  cout << "Unknown dzpuis in sol_poisson_tau ..." << endl ;
370  abort() ;
371  }
372 
373  ligne_courant += nbr_cl ;
374 
375  // Multiplication of the source :
376  double indic = 1 ;
377  switch (dzpuis) {
378  case 4 :
379  indic = alpha*alpha ;
380  break ;
381  case 3 :
382  indic = alpha ;
383  break ;
384  default :
385  break ;
386  }
387 
388  // L'operateur :
389  for (int lig=0 ; lig<nr-1-nbr_cl ; lig++) {
390  for (int col=0 ; col<nr ; col++)
391  systeme.set(lig+ligne_courant,col+column_courant) = (*work)(lig,col) ;
392  sec_membre.set(lig+ligne_courant) = indic*source(nz-1, k, j, lig) ;
393  }
394  delete work ;
395 
396  // Solving the system:
397  systeme.set_band (max_nr, max_nr) ;
398  systeme.set_lu() ;
399  Tbl soluce (systeme.inverse(sec_membre)) ;
400 
401  // On range :
402  int conte = 0 ;
403  for (int l=0 ; l<nz ; l++) {
404  nr = mapping.get_mg()->get_nr(l) ;
405  for (int i=0 ; i<nr ; i++) {
406  resultat.set(l,k,j,i) = soluce(conte) ;
407  conte ++ ;
408  }
409  }
410 
411  }
412 
413  return resultat ;
414 }
415 }
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
Definition: mtbl_cf.h:453
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
Lorene prototypes.
Definition: app_hor.h:64