LORENE
vector_divfree_A_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 vector_divfree_A_tau[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/vector_divfree_A_tau.C,v 1.4 2014/10/13 08:53:45 j_novak Exp $" ;
24 
25 /*
26  * $Id: vector_divfree_A_tau.C,v 1.4 2014/10/13 08:53:45 j_novak Exp $
27  * $Log: vector_divfree_A_tau.C,v $
28  * Revision 1.4 2014/10/13 08:53:45 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.3 2014/10/06 15:13:20 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.2 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.1 2008/08/27 09:01:27 jl_cornou
39  * Methods for solving Dirac systems for divergence free vectors
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/Tensor/vector_divfree_A_tau.C,v 1.4 2014/10/13 08:53:45 j_novak Exp $
56  *
57  */
58 
59 // C headers
60 #include <cassert>
61 #include <cmath>
62 
63 // Lorene headers
64 #include "tensor.h"
65 #include "diff.h"
66 #include "proto.h"
67 #include "param.h"
68 #include "matrice.h"
69 #include "type_parite.h"
70 
71 
72 namespace Lorene {
73 void Vector_divfree::sol_Dirac_A_tau(const Scalar& aaa, Scalar& eta_tilde, Scalar& vr, const Param* par_bc) const {
74 
75  const Map_af* mp_aff = dynamic_cast<const Map_af*>(mp) ;
76  assert(mp_aff != 0x0) ; // Only affine mapping for the moment
77 
78  const Mg3d& mgrid = *mp_aff->get_mg();
79  assert((mgrid.get_type_r(0) == RARE) || (mgrid.get_type_r(0) == FIN));
80  if (aaa.get_etat() == ETATZERO) {
81  eta_tilde = 0;
82  vr = 0 ;
83  return ;
84  }
85  assert(aaa.get_etat() != ETATNONDEF) ;
86  int nz = mgrid.get_nzone();
87  int nzm1 = nz - 1 ;
88  bool ced = (mgrid.get_type_r(nzm1) == UNSURR) ;
89  int n_shell = ced ? nz-2 : nzm1 ;
90  int nz_bc = nzm1 ;
91  if (par_bc != 0x0)
92  if (par_bc->get_n_int() > 0) nz_bc = par_bc->get_int();
93  n_shell = (nz_bc < n_shell ? nz_bc : n_shell) ;
94  //bool cedbc (ced && (nz_bc == nzm1));
95 
96 //#ifndef NDEBUG
97 // if (!cedbc) {
98 // assert(par_bc != 0x0) ;
99 // assert(par_bc->get_n_tbl_mod() >= 3) ;
100 // }
101 //#endif
102  int nt = mgrid.get_nt(0) ;
103  int np = mgrid.get_np(0) ;
104  int nr ;
105 
106  Scalar source = aaa ;
107  Scalar source_coq = aaa ;
108  source_coq.annule_domain(0);
109 
110  int dzpuis = source.get_dzpuis();
111 
112  if (ced) source_coq.annule_domain(nzm1) ;
113  source_coq.mult_r() ;
114  source.set_spectral_va().ylm() ;
115  source_coq.set_spectral_va().ylm() ;
116  Base_val base = source.get_spectral_base() ;
117  base.mult_x() ;
118 
119  eta_tilde.annule_hard() ;
120  eta_tilde.set_spectral_base(base) ;
121  eta_tilde.set_spectral_va().set_etat_cf_qcq() ;
122  eta_tilde.set_spectral_va().c_cf->annule_hard() ;
123  vr.annule_hard() ;
124  vr.set_spectral_base(base) ;
125  vr.set_spectral_va().set_etat_cf_qcq() ;
126  vr.set_spectral_va().c_cf->annule_hard() ;
127  Mtbl_cf& meta = *eta_tilde.set_spectral_va().c_cf ;
128  Mtbl_cf& mvr = *vr.set_spectral_va().c_cf ;
129 
130 
131  int base_r ;
132 
133 
134  // Determination of the size of the systeme :
135  int size = 0 ;
136  int max_nr = 0 ;
137  for (int l=0 ; l<nz ; l++) {
138  nr = mgrid.get_nr(l) ;
139  size += 2*nr ;
140  if (nr > max_nr)
141  max_nr = nr ;
142  }
143 
144  Matrice systeme (size, size) ;
145  systeme.set_etat_qcq() ;
146  Tbl sec_membre (size) ;
147 
148  np = mgrid.get_np(0) ;
149  nt = mgrid.get_nt(0) ;
150 
151 
152  double alpha, beta ;
153  int l_quant, m_quant ;
154 
155  // On bosse pour chaque l, m :
156  for (int k=0 ; k<np+1 ; k++)
157  for (int j=0 ; j<nt ; j++)
158  if (nullite_plm(j, nt, k, np, base) == 1) {
159 
160  systeme.annule_hard() ;
161  sec_membre.annule_hard() ;
162 
163  int column_courant = 0 ;
164  int ligne_courant = 0 ;
165 
166  //--------------------------
167  // NUCLEUS
168  //--------------------------
169 
170  nr = mgrid.get_nr(0) ;
171  alpha = (*mp_aff).get_alpha()[0] ;
172  base.give_quant_numbers (0, k, j, m_quant, l_quant, base_r) ;
173  Diff_dsdx odn(base_r, nr) ; const Matrice& mdn = odn.get_matrice() ;
174  Diff_sx osn(base_r, nr) ; const Matrice& msn = osn.get_matrice() ;
175 
176 
177 
178  int nbr_cl = 0 ;
179  // RARE case
180  if (source.get_mp().get_mg()->get_type_r(0) == RARE) {
181  // regularity conditions for eta :
182  if (l_quant > 1) {
183  nbr_cl = 1 ;
184  if (l_quant%2==0) {
185  //Even case
186  for (int col=0 ; col<nr ; col++) {
187  if (col%2==0) {
188  systeme.set(ligne_courant, col+column_courant) = 1 ;
189  systeme.set(ligne_courant, col+column_courant+nr) = 0 ; }
190  else {
191  systeme.set(ligne_courant, col+column_courant) = -1 ;
192  systeme.set(ligne_courant, col+column_courant+nr) = 0 ; }
193  }
194  }
195  else {
196  //Odd case
197  for (int col=0 ; col<nr ; col++) {
198  if (col%2==0) {
199  systeme.set(ligne_courant, col+column_courant) = 2*col+1 ;
200  systeme.set(ligne_courant, col+column_courant+nr) = 0 ; }
201  else {
202  systeme.set(ligne_courant, col+column_courant) = -(2*col+1) ;
203  systeme.set(ligne_courant, col+column_courant+nr) = 0 ; }
204  }
205  }
206  }
207  }
208 
209  // R_JACO02 case
210  else {
211  assert (base_r == R_JACO02) ;
212  // regularity conditions for eta :
213  if (l_quant == 0) {
214  nbr_cl = 1 ;
215  for (int col=0 ; col<nr ; col++) {
216  systeme.set(ligne_courant, col+column_courant) = col*(col+1)*(col+2)*(col+3)/double(12)*(2*(col%2)-1);
217  systeme.set(ligne_courant, col+column_courant+nr) = 0 ;
218  }
219  }
220  else if (l_quant == 1) {
221  nbr_cl = 1 ;
222  for (int col=0 ; col<nr ; col++) {
223  systeme.set(ligne_courant, col+column_courant) = (col+1)*(col+2)/double(2)*(1-2*(col%2)) ;
224  systeme.set(ligne_courant, col+column_courant+nr) = 0 ;
225  }
226  }
227  else {
228  nbr_cl = 2 ;
229  for (int col=0 ; col<nr ; col++) {
230  systeme.set(ligne_courant, col+column_courant) = (col+1)*(col+2)/double(2)*(1-2*(col%2)) ;
231  systeme.set(ligne_courant, col+column_courant+nr) = 0 ;
232  systeme.set(ligne_courant+1, col+column_courant) = col*(col+1)*(col+2)*(col+3)/double(12)*(2*(col%2)-1) ;
233  systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ;
234  }
235  }
236  }
237  ligne_courant += nbr_cl ;
238 
239  // Premiere partie de l'operateur :
240  for (int lig=0 ; lig<nr-1-nbr_cl ; lig++) {
241  for (int col=0 ; col<nr ; col++) {
242  systeme.set(lig+ligne_courant,col+column_courant) = mdn(lig,col) + msn(lig,col) ;
243  sec_membre.set(lig+ligne_courant) = alpha*alpha*(*source.get_spectral_va().c_cf)(0, k, j, lig) ;
244  systeme.set(lig+ligne_courant,col+column_courant+nr)= - msn(lig,col) ;
245  // systeme.set(lig+ligne_courant+nr) = 0 ;
246  }
247  }
248 
249 
250  ligne_courant += nr-1-nbr_cl ;
251 
252 
253  // RARE case
254  if (source.get_mp().get_mg()->get_type_r(0) == RARE) {
255  // regularity conditions for vr :
256  if (l_quant > 1) {
257  nbr_cl = 1 ;
258  if (l_quant%2==0) {
259  //Even case
260  for (int col=0 ; col<nr ; col++) {
261  if (col%2==0) {
262  systeme.set(ligne_courant, col+column_courant) = 0 ;
263  systeme.set(ligne_courant, col+column_courant+nr) = 1 ; }
264  else {
265  systeme.set(ligne_courant, col+column_courant) = 0 ;
266  systeme.set(ligne_courant, col+column_courant+nr) = -1 ; }
267  }
268  }
269 
270  else {
271  //Odd case
272  for (int col=0 ; col<nr ; col++) {
273  if (col%2==0) {
274  systeme.set(ligne_courant, col+column_courant) = 0 ;
275  systeme.set(ligne_courant, col+column_courant+nr) = 2*col+1 ; }
276  else {
277  systeme.set(ligne_courant, col+column_courant) = 0 ;
278  systeme.set(ligne_courant, col+column_courant+nr) = -(2*col+1) ; }
279  }
280  }
281  }
282  }
283 
284  // R_JACO02 case
285  else {
286  assert (base_r == R_JACO02) ;
287  // regularity conditions for vr :
288  if (l_quant == 0) {
289  nbr_cl = 1 ;
290  for (int col=0 ; col<nr ; col++) {
291  systeme.set(ligne_courant, col+column_courant) = 0 ;
292  systeme.set(ligne_courant, col+column_courant+nr) = col*(col+1)*(col+2)*(col+3)/double(12)*(2*(col%2)-1) ;
293  }
294  }
295  else if (l_quant == 1) {
296  nbr_cl = 1 ;
297  for (int col=0 ; col<nr ; col++) {
298  systeme.set(ligne_courant, col+column_courant) = 0 ;
299  systeme.set(ligne_courant, col+column_courant+nr) = (col+1)*(col+2)/double(2)*(1-2*(col%2)) ;
300  }
301  }
302  else {
303  nbr_cl = 2 ;
304  for (int col=0 ; col<nr ; col++) {
305  systeme.set(ligne_courant, col+column_courant) = 0 ;
306  systeme.set(ligne_courant, col+column_courant+nr) = (col+1)*(col+2)/double(2)*(1-2*(col%2)) ;
307  systeme.set(ligne_courant+1, col+column_courant) = 0 ;
308  systeme.set(ligne_courant+1, col+column_courant+nr) = col*(col+1)*(col+2)*(col+3)/double(12)*(2*(col%2)-1) ;
309  }
310  }
311  }
312  ligne_courant += nbr_cl ;
313 
314  // Deuxieme partie de l'operateur
315  for (int lig=0 ; lig<nr-1-nbr_cl ; lig++) {
316  for (int col=0 ; col<nr ; col++) {
317  systeme.set(lig+ligne_courant,col+column_courant) = - l_quant*(l_quant+1)*msn(lig,col) ;
318  sec_membre.set(lig+ligne_courant) = 0 ;
319  systeme.set(lig+ligne_courant,col+column_courant+nr)= mdn(lig,col) + 2*msn(lig,col) ;
320  // systeme.set(lig+ligne_courant+nr) = 0 ;
321  }
322  }
323 
324 
325  ligne_courant += nr-1-nbr_cl ;
326 
327 
328  // Le raccord pour eta
329  for (int col=0 ; col<nr ; col++) {
330  if (source.get_mp().get_mg()->get_type_r(0) == RARE) {
331  // La fonction eta
332  systeme.set(ligne_courant, col+column_courant) = 1 ;
333  systeme.set(ligne_courant, col+column_courant+nr) = 0 ;
334  // Sa derivee :
335  if (l_quant%2==0) {
336  systeme.set(ligne_courant+1, col+column_courant) = 4*col*col/alpha ;
337  systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ;
338  }
339  else {
340  systeme.set(ligne_courant+1, col+column_courant) = (2*col+1)*(2*col+1)/alpha ;
341  systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ;
342  }
343  }
344  else { // Cas R_JACO02
345  assert( base_r == R_JACO02 ) ;
346  // La fonction eta
347  systeme.set(ligne_courant, col+column_courant) = 1 ;
348  systeme.set(ligne_courant, col+column_courant+nr) = 0 ;
349  // Sa derivee :
350  systeme.set(ligne_courant+1, col+column_courant) = col*(col+3)/double(2)/alpha ;
351  systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ;
352  }
353  }
354 
355  ligne_courant += 2 ;
356  // Le raccord pour vr
357  for (int col=0 ; col<nr ; col++) {
358  if (source.get_mp().get_mg()->get_type_r(0) == RARE) {
359  // La fonction vr
360  systeme.set(ligne_courant, col+column_courant) = 0 ;
361  systeme.set(ligne_courant, col+column_courant+nr) = 1 ;
362  // Sa derivee :
363  if (l_quant%2==0) {
364  systeme.set(ligne_courant+1, col+column_courant) = 0 ;
365  systeme.set(ligne_courant+1, col+column_courant+nr) = 4*col*col/alpha ;
366  }
367  else {
368  systeme.set(ligne_courant+1, col+column_courant) = 0 ;
369  systeme.set(ligne_courant+1, col+column_courant+nr) = (2*col+1)*(2*col+1)/alpha ;
370  }
371  }
372  else { // Cas R_JACO02
373  assert( base_r == R_JACO02 ) ;
374  // La fonction vr
375  systeme.set(ligne_courant, col+column_courant) = 0 ;
376  systeme.set(ligne_courant, col+column_courant+nr) = 1 ;
377  // Sa derivee :
378  systeme.set(ligne_courant+1, col+column_courant) = 0 ;
379  systeme.set(ligne_courant+1, col+column_courant+nr) = col*(col+3)/double(2)/alpha ;
380  }
381  }
382 
383  ligne_courant -= 2 ; // On retourne pour le raccord dans la prochaine zone
384 
385  column_courant += 2*nr ; // On va changer de zone
386 
387  //--------------------------
388  // SHELLS
389  //--------------------------
390  for (int l=1 ; l<nz-1 ; l++) {
391 
392  nr = mgrid.get_nr(l) ;
393  alpha = (*mp_aff).get_alpha()[l] ;
394  beta = (*mp_aff).get_beta()[l] ;
395 
396  base.give_quant_numbers (l, k, j, m_quant, l_quant, base_r) ;
397  Diff_id ois(base_r, nr) ; const Matrice& mis = ois.get_matrice() ;
398  Diff_xdsdx oxds(base_r, nr) ; const Matrice& mxds = oxds.get_matrice() ;
399 
400  // matching with previous domain :
401  for (int col=0 ; col<nr ; col++) {
402  // La fonction eta
403  if (col%2==0) {
404  systeme.set(ligne_courant, col+column_courant) = -1 ;
405  systeme.set(ligne_courant, col+column_courant+nr) = 0 ; }
406  else {
407  systeme.set(ligne_courant, col+column_courant) = 1 ;
408  systeme.set(ligne_courant, col+column_courant+nr) = 0 ; }
409  // Sa derivee :
410  if (col%2==0) {
411  systeme.set(ligne_courant+1, col+column_courant) = col*col/alpha ;
412  systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ; }
413  else {
414  systeme.set(ligne_courant+1, col+column_courant) = -col*col/alpha ;
415  systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ; }
416  }
417  ligne_courant += 2 ;
418 
419  // matching with previous domain :
420  for (int col=0 ; col<nr ; col++) {
421  // La fonction vr
422  if (col%2==0) {
423  systeme.set(ligne_courant, col+column_courant) = 0 ;
424  systeme.set(ligne_courant, col+column_courant+nr) = -1 ; }
425  else {
426  systeme.set(ligne_courant, col+column_courant) = 0 ;
427  systeme.set(ligne_courant, col+column_courant+nr) = 1 ; }
428  // Sa derivee :
429  if (col%2==0) {
430  systeme.set(ligne_courant+1, col+column_courant) = 0 ;
431  systeme.set(ligne_courant+1, col+column_courant+nr) = col*col/alpha ; }
432  else {
433  systeme.set(ligne_courant+1, col+column_courant) = 0 ;
434  systeme.set(ligne_courant+1, col+column_courant+nr) = -col*col/alpha ; }
435  }
436  ligne_courant += 2 ;
437 
438 
439  // L'operateur (premiere et deuxieme parties) :
440 
441  // source must be multiplied by (x+echelle)^2
442  Tbl source_aux(nr) ;
443  source_aux.set_etat_qcq() ;
444  for (int i=0 ; i<nr ; i++)
445  source_aux.set(i) = (*source.get_spectral_va().c_cf)(l,k,j,i)*alpha*alpha ;
446  Tbl xso(source_aux) ;
447  Tbl xxso(source_aux) ;
448  multx_1d(nr, &xso.t, R_CHEB) ;
449  multx_1d(nr, &xxso.t, R_CHEB) ;
450  multx_1d(nr, &xxso.t, R_CHEB) ;
451  source_aux = beta*beta/alpha/alpha*source_aux+2*beta/alpha*xso+xxso ;
452 
453  for (int lig=0 ; lig<nr-2 ; lig++) {
454  for (int col=0 ; col<nr ; col++) {
455  systeme.set(lig+ligne_courant,col+column_courant) = mxds(lig,col) + mis(lig,col) ;
456  systeme.set(lig+ligne_courant,col+column_courant+nr) = -mis(lig,col) ;
457  sec_membre.set(lig+ligne_courant) = source_aux(lig) ;
458  systeme.set(lig+ligne_courant+nr-2, col+column_courant) = -l_quant*(l_quant+1) ;
459  systeme.set(lig+ligne_courant+nr-2, col+column_courant+nr) = mxds(lig,col) + 2*mis(lig,col) ;
460  sec_membre.set(lig+ligne_courant+nr-2) = 0 ; }
461  }
462 
463 
464  ligne_courant += 2*nr-4 ;
465  // Matching with the next domain :
466  for (int col=0 ; col<nr ; col++) {
467  // La fonction eta
468  systeme.set(ligne_courant, col+column_courant) = 1 ;
469  systeme.set(ligne_courant, col+column_courant+nr) = 0 ;
470  // Sa derivee :
471  systeme.set(ligne_courant+1, col+column_courant) = col*col/alpha ;
472  systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ;
473  // La fonction vr
474  systeme.set(ligne_courant+2, col+column_courant) = 0 ;
475  systeme.set(ligne_courant+2, col+column_courant+nr) = 1 ;
476  // Sa derivee
477  systeme.set(ligne_courant+3, col+column_courant) = 0 ;
478  systeme.set(ligne_courant+3, col+column_courant+nr) = col*col/alpha ;
479  }
480 
481  column_courant += 2*nr ;
482  }
483 
484 
485  //--------------------------
486  // ZEC
487  //--------------------------
488  nr = mgrid.get_nr(nz-1) ;
489  alpha = (*mp_aff).get_alpha()[nz-1] ;
490  beta = (*mp_aff).get_beta()[nz-1] ;
491 
492  base.give_quant_numbers (nz-1, k, j, m_quant, l_quant, base_r) ;
493  //work = new Matrice(laplacien_mat(nr, l_quant, 0., dzpuis, base_r)) ;
494  Diff_dsdx odzec(base_r, nr) ; const Matrice& mdzec = odzec.get_matrice() ;
495  Diff_sx oszec(base_r, nr) ; const Matrice& mszec = oszec.get_matrice() ;
496 
497 
498  // Matching with the previous domain :
499  for (int col=0 ; col<nr ; col++) {
500  // La fonction eta
501  if (col%2==0) {
502  systeme.set(ligne_courant, col+column_courant) = -1 ;
503  systeme.set(ligne_courant, col+column_courant+nr) = 0 ; }
504  else {
505  systeme.set(ligne_courant, col+column_courant) = 1 ;
506  systeme.set(ligne_courant, col+column_courant+nr) = 0 ; }
507  // Sa derivee :
508  if (col%2==0) {
509  systeme.set(ligne_courant+1, col+column_courant) = -4*alpha*col*col ;
510  systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ; }
511  else {
512  systeme.set(ligne_courant+1, col+column_courant) = 4*alpha*col*col ;
513  systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ; }
514  // La fonction vr
515  if (col%2==0) {
516  systeme.set(ligne_courant+2, col+column_courant) = -1 ;
517  systeme.set(ligne_courant+2, col+column_courant+nr) = 0 ; }
518  else {
519  systeme.set(ligne_courant+2, col+column_courant) = 1 ;
520  systeme.set(ligne_courant+2, col+column_courant+nr) = 0 ; }
521  // Sa derivee :
522  if (col%2==0) {
523  systeme.set(ligne_courant+3, col+column_courant) = -4*alpha*col*col ;
524  systeme.set(ligne_courant+3, col+column_courant+nr) = 0 ; }
525  else {
526  systeme.set(ligne_courant+3, col+column_courant) = 4*alpha*col*col ;
527  systeme.set(ligne_courant+3, col+column_courant+nr) = 0 ; }
528  }
529  ligne_courant += 4 ;
530 
531  // Regularity and BC at infinity ?
532  nbr_cl =0 ;
533  switch (dzpuis) {
534  case 4 :
535  if (l_quant==0) {
536  nbr_cl = 1 ;
537  // Only BC at infinity :
538  for (int col=0 ; col<nr ; col++) {
539  systeme.set(ligne_courant, col+column_courant) = 1 ;
540  systeme.set(ligne_courant, col+column_courant+nr) = 0 ;
541  systeme.set(ligne_courant+1, col+column_courant) = 0 ;
542  systeme.set(ligne_courant+1, col+column_courant+nr) = 1 ; }
543  }
544  else {
545  nbr_cl = 2 ;
546  // BC at infinity :
547  for (int col=0 ; col<nr ; col++) {
548  systeme.set(ligne_courant, col+column_courant) = 1 ;
549  systeme.set(ligne_courant, col+column_courant+nr) = 0 ;
550  systeme.set(ligne_courant+1, col+column_courant) = 0 ;
551  systeme.set(ligne_courant+1, col+column_courant+nr) = 1; }
552  // Regularity :
553  for (int col=0 ; col<nr ; col++) {
554  systeme.set(ligne_courant+2, col+column_courant) = -4*alpha*col*col ;
555  systeme.set(ligne_courant+2, col+column_courant+nr) = 0 ;
556  systeme.set(ligne_courant+3, col+column_courant) = 0 ;
557  systeme.set(ligne_courant+3, col+column_courant+nr) = -4*alpha*col*col ; }
558  }
559  break ;
560 
561  case 3 :
562  nbr_cl = 1 ;
563  // Only BC at infinity :
564  for (int col=0 ; col<nr ; col++) {
565  systeme.set(ligne_courant, col+column_courant) = 1 ;
566  systeme.set(ligne_courant, col+column_courant+nr) = 0 ;
567  systeme.set(ligne_courant+1, col+column_courant) = 0 ;
568  systeme.set(ligne_courant+1, col+column_courant+nr) = 1 ; }
569  break ;
570 
571  case 2 :
572  if (l_quant==0) {
573  nbr_cl = 1 ;
574  // Only BC at infinity :
575  for (int col=0 ; col<nr ; col++) {
576  systeme.set(ligne_courant, col+column_courant) = 1 ;
577  systeme.set(ligne_courant, col+column_courant+nr) = 0 ;
578  systeme.set(ligne_courant+1, col+column_courant) = 0 ;
579  systeme.set(ligne_courant+1, col+column_courant+1) = 1 ; }
580  }
581  break ;
582  default :
583  cout << "Unknown dzpuis in vector_divfree_A ..." << endl ;
584  abort() ;
585  }
586 
587  ligne_courant += nbr_cl ;
588 
589  // Multiplication of the source :
590  double indic = 1 ;
591  switch (dzpuis) {
592  case 4 :
593  indic = alpha*alpha ;
594  break ;
595  case 3 :
596  indic = alpha ;
597  break ;
598  default :
599  break ;
600  }
601 
602  // L'operateur :
603  for (int lig=0 ; lig<nr-1-nbr_cl ; lig++) {
604  for (int col=0 ; col<nr ; col++) {
605  systeme.set(lig+ligne_courant,col+column_courant) = mdzec(lig,col)+mszec(lig,col) ;
606  systeme.set(lig+ligne_courant,col+column_courant+nr) = -mszec(lig,col) ;
607  sec_membre.set(lig+ligne_courant) = indic*(*source.get_spectral_va().c_cf)(nz-1, k, j, lig) ;
608  systeme.set(lig+ligne_courant+nr-1-nbr_cl,col+column_courant)= - l_quant*(l_quant+1)*mszec(lig,col) ;
609  systeme.set(lig+ligne_courant+nr-1-nbr_cl,col+column_courant+nr) = mdzec(lig,col)+2*mszec(lig,col) ;
610  sec_membre.set(lig+ligne_courant+nr-1-nbr_cl) = 0 ; }
611  }
612 
613 
614  // Solving the system:
615  systeme.set_band (max_nr, max_nr) ;
616  systeme.set_lu() ;
617  Tbl soluce (systeme.inverse(sec_membre)) ;
618 
619  // On range :
620  int conte = 0 ;
621  for (int l=0 ; l<nz ; l++) {
622  nr = mgrid.get_nr(l) ;
623  for (int i=0 ; i<nr ; i++) {
624  meta.set(l,k,j,i) = soluce(conte) ;
625  mvr.set(l,k,i,j) = soluce(conte+nr) ;
626  conte ++ ;
627  }
628  }
629 }
630 if (eta_tilde.set_spectral_va().c != 0x0)
631  delete eta_tilde.set_spectral_va().c ;
632  eta_tilde.set_spectral_va().c = 0x0 ;
633  eta_tilde.set_spectral_va().ylm_i() ;
634 
635  if (vr.set_spectral_va().c != 0x0)
636  delete vr.set_spectral_va().c ;
637  vr.set_spectral_va().c = 0x0 ;
638  vr.set_spectral_va().ylm_i() ;
639 }
640 }
Lorene::Scalar::get_spectral_base
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Definition: scalar.h:1328
Lorene::Diff_sx::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_sx.C:103
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
Lorene::Scalar::mult_r
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
Definition: scalar_r_manip.C:211
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:134
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
Lorene::Valeur::set_etat_cf_qcq
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
Definition: valeur.C:715
Lorene::Base_val::give_quant_numbers
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Definition: base_val_quantum.C:84
Lorene::Matrice::set_lu
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition: matrice.C:395
Lorene
Lorene prototypes.
Definition: app_hor.h:67
Lorene::Mtbl_cf::set
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:304
Lorene::Matrice::inverse
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition: matrice.C:427
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:141
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:777
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Lorene::Diff_dsdx::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_dsdx.C:97
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Lorene::Diff_dsdx
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:129
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
R_JACO02
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
Lorene::Scalar::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:386
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
Lorene::Diff_id
Class for the elementary differential operator Identity (see the base class Diff ).
Definition: diff.h:210
Lorene::Param::get_n_int
int get_n_int() const
Returns the number of stored int &#39;s addresses.
Definition: param.C:242
Lorene::Param::get_int
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:295
Lorene::Matrice
Matrix handling.
Definition: matrice.h:152
Lorene::Scalar::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:563
Lorene::Tbl::t
double * t
The array of double.
Definition: tbl.h:176
Lorene::Valeur::c
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:309
Lorene::Param
Parameter storage.
Definition: param.h:125
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Lorene::Matrice::annule_hard
void annule_hard()
Sets the logical state to ETATQCQ (undefined state).
Definition: matrice.C:196
Lorene::Matrice::set
double & set(int j, int i)
Read/write of a particuliar element.
Definition: matrice.h:277
Lorene::Matrice::set_band
void set_band(int up, int low) const
Calculate the band storage of *std.
Definition: matrice.C:367
Lorene::Scalar::set_spectral_base
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
Definition: scalar.C:803
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Lorene::Mg3d
Multi-domain grid.
Definition: grilles.h:279
Lorene::Base_val
Bases of the spectral expansions.
Definition: base_val.h:325
Lorene::Mtbl_cf::annule_hard
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition: mtbl_cf.C:315
Lorene::Map_af
Affine radial mapping.
Definition: map.h:2042
Lorene::Diff_xdsdx
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:409
Lorene::Base_val::mult_x
void mult_x()
The basis is transformed as with a multiplication by .
Definition: base_val_manip.C:160
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:196
Lorene::Matrice::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:178
Lorene::Tbl
Basic array class.
Definition: tbl.h:164
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
Lorene::Diff_id::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_id.C:94
Lorene::Diff_sx
Class for the elementary differential operator division by (see the base class Diff )...
Definition: diff.h:329
Lorene::Mg3d::get_type_r
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:491
Lorene::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:301
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:874
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:375
Lorene::Diff_xdsdx::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_xdsdx.C:101
Lorene::Scalar::get_spectral_va
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607
Lorene::Vector_divfree::sol_Dirac_A_tau
void sol_Dirac_A_tau(const Scalar &aaa, Scalar &eta, Scalar &vr, const Param *par_bc=0x0) const
Solves via a tau method a system of two-coupled first-order PDEs obtained from the divergence-free co...
Definition: vector_divfree_A_tau.C:73
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166