LORENE
map_af_poisson_regu.C
1 /*
2  * Method of the class Map_af for the resolution of the scalar Poisson
3  * equation by using regularized source.
4  */
5 
6 /*
7  * Copyright (c) 2000-2001 Keisuke Taniguchi
8  *
9  * This file is part of LORENE.
10  *
11  * LORENE is free software; you can redistribute it and/or modify
12  * it under the terms of the GNU General Public License as published by
13  * the Free Software Foundation; either version 2 of the License, or
14  * (at your option) any later version.
15  *
16  * LORENE is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  * GNU General Public License for more details.
20  *
21  * You should have received a copy of the GNU General Public License
22  * along with LORENE; if not, write to the Free Software
23  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24  *
25  */
26 
27 
28 char map_af_poisson_regu_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_af_poisson_regu.C,v 1.5 2014/10/13 08:53:03 j_novak Exp $" ;
29 
30 /*
31  * $Id: map_af_poisson_regu.C,v 1.5 2014/10/13 08:53:03 j_novak Exp $
32  * $Log: map_af_poisson_regu.C,v $
33  * Revision 1.5 2014/10/13 08:53:03 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.4 2014/10/06 15:13:12 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.3 2003/12/19 16:21:43 j_novak
40  * Shadow hunt
41  *
42  * Revision 1.2 2003/10/03 15:58:48 j_novak
43  * Cleaning of some headers
44  *
45  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
46  * LORENE
47  *
48  * Revision 2.13 2000/10/25 16:05:11 keisuke
49  * Remake for the arbitrary regularization degree (k_div).
50  *
51  * Revision 2.12 2000/10/06 15:31:41 keisuke
52  * Suppress assertions for nilm_cos and nilm_sin.
53  * Add warnings for nilm_cos and nilm_sin.
54  *
55  * Revision 2.11 2000/09/25 15:02:48 keisuke
56  * modify the derivative duu_div.
57  *
58  * Revision 2.10 2000/09/11 13:50:57 keisuke
59  * Set the basis of duu_div to the spherical one.
60  *
61  * Revision 2.9 2000/09/11 10:15:06 keisuke
62  * Change the method to reconstruct source_regu.
63  *
64  * Revision 2.8 2000/09/09 14:51:20 keisuke
65  * Suppress uu_regu.set_dzpuis(0).
66  *
67  * Revision 2.7 2000/09/07 15:29:50 keisuke
68  * Add a new argument Cmp& uu.
69  *
70  * Revision 2.6 2000/09/06 10:28:42 keisuke
71  * Modify the scaling for derivatives.
72  *
73  * Revision 2.5 2000/09/04 15:55:38 keisuke
74  * Include the polar and azimuthal parts of duu_div.
75  *
76  * Revision 2.4 2000/09/04 13:08:39 keisuke
77  * Change alpha[0] into mp_radial->dxdr.
78  *
79  * Revision 2.3 2000/09/01 08:55:55 keisuke
80  * Change val_r into alpha[0].
81  *
82  * Revision 2.2 2000/08/31 15:59:51 keisuke
83  * Modify the arguments.
84  *
85  * Revision 2.1 2000/08/28 16:11:44 keisuke
86  * Add "int nzet" in the argumant.
87  *
88  * Revision 2.0 2000/08/25 08:48:36 keisuke
89  * *** empty log message ***
90  *
91  *
92  * $Header: /cvsroot/Lorene/C++/Source/Map/map_af_poisson_regu.C,v 1.5 2014/10/13 08:53:03 j_novak Exp $
93  *
94  */
95 
96 // headers C
97 #include <cmath>
98 
99 // Header Lorene
100 #include "tenseur.h"
101 #include "matrice.h"
102 #include "param.h"
103 #include "proto.h"
104 
105 //******************************************************************
106 
107 namespace Lorene {
108 
109 void Map_af::poisson_regular(const Cmp& source, int k_div, int nzet,
110  double unsgam1, Param& par, Cmp& uu,
111  Cmp& uu_regu, Cmp& uu_div, Tenseur& duu_div,
112  Cmp& source_regu, Cmp& source_div) const {
113 
114 
115  assert(source.get_etat() != ETATNONDEF) ;
116  assert(source.get_mp()->get_mg() == mg) ;
117  assert(k_div > 0) ;
118 
119  double aa = unsgam1 ; // exponent of the specific enthalpy
120 
121  int nzm1 = mg->get_nzone() - 1;
122  int nr = mg->get_nr(0) ;
123  int nt = mg->get_nt(0) ;
124  int np = mg->get_np(0) ;
125 
126  assert(nr-k_div > 0) ;
127 
128  // --------------------------------------------
129  // Expansion of "source" by T_i Y_l^m
130  // --------------------------------------------
131 
132  // Expansion by cos(j \theta) and sin(j \theta) for theta direction
133  // ----------------------------------------------------------------
134 
135  const Valeur& sourva = source.va ;
136  assert(sourva.get_etat() == ETATQCQ) ;
137 
138  Valeur rho(sourva.get_mg()) ;
139  sourva.coef() ;
140  rho = *(sourva.c_cf) ;
141 
142  // Expansion by Legendre function for theta direction
143  // --------------------------------------------------
144 
145  rho.ylm() ;
146 
147  // Obtaining the coefficients of the given source
148  // ----------------------------------------------
149 
150  Tbl& ccf = *((rho.c_cf)->t[0]) ;
151 
152  Tbl nilm_cos(np/2+1, 2*nt, nr) ;
153  Tbl nilm_sin(np/2+1, 2*nt, nr) ;
154 
155  nilm_cos.set_etat_qcq() ;
156  nilm_sin.set_etat_qcq() ;
157 
158  for (int k=0; k<=np; k+=2) {
159 
160  int m = k / 2 ;
161 
162  for (int j=0; j<nt; j++) {
163 
164  int l ;
165  if(m%2 == 0) {
166  l = 2 * j ;
167  }
168  else
169  l = 2 * j + 1 ;
170 
171  for (int i=0; i<nr; i++) {
172 
173  nilm_cos.set(m, l, i) = ccf(k, j, i) ;
174  nilm_sin.set(m, l, i) = ccf(k+1, j, i) ;
175 
176  }
177  }
178  }
179 
180  // -----------------------------------------------------
181  // Expansion of "analytic source" by T_i Y_l^m
182  // -----------------------------------------------------
183 
184  const Grille3d& grid = *(mg->get_grille3d(0)) ;
185 
186  Tbl cf_cil(2*nt, nr, k_div) ;
187  cf_cil.set_etat_qcq() ;
188 
189  int deg[3] ;
190  int dim[3] ;
191 
192  deg[0] = 1 ;
193  deg[1] = 1 ;
194  deg[2] = nr ;
195  dim[0] = 1 ;
196  dim[1] = 1 ;
197  dim[2] = nr ;
198 
199  double* tmp1 = new double[nr] ;
200 
201  for (int k_deg=1; k_deg<=k_div; k_deg++) {
202 
203  for (int l=0; l<2*nt; l++) {
204 
205  for (int i=0; i<nr; i++) {
206 
207  double xi = grid.x[i] ;
208 
209  tmp1[i] = (aa + k_deg + 1.) *
210  ( -(4. * l + 6.) * pow(1. - xi * xi, aa + k_deg) * pow(xi, l)
211  + 4. * (aa + k_deg) * pow(1. - xi * xi, aa + k_deg - 1.) *
212  pow(xi, l + 2.) ) ;
213 
214  }
215 
216  if (l%2 == 0) {
217  cfrchebp(deg, dim, tmp1, dim, tmp1) ;
218  }
219  else
220  cfrchebi(deg, dim, tmp1, dim, tmp1) ;
221 
222  for (int i=0; i<nr; i++) {
223 
224  cf_cil.set(l, i, k_deg-1) = tmp1[i] ;
225 
226  }
227  }
228  }
229 
230  // Calculation of the coefficients : Solve the simultaneous equations
231  // -------------------------------
232 
233  Tbl alm_cos(np/2+1, 2*nt, k_div) ;
234  Tbl alm_sin(np/2+1, 2*nt, k_div) ;
235 
236  alm_cos.set_etat_qcq() ;
237  alm_sin.set_etat_qcq() ;
238 
239  Matrice matrix(k_div, k_div) ;
240  matrix.set_etat_qcq() ;
241 
242  Tbl rhs_cos(k_div) ;
243  Tbl rhs_sin(k_div) ;
244 
245  rhs_cos.set_etat_qcq() ;
246  rhs_sin.set_etat_qcq() ;
247 
248  for (int k=0; k<=np; k+=2) {
249 
250  int m = k / 2 ;
251 
252  for (int j=0; j<nt; j++) {
253 
254  int l ;
255  if(m%2 == 0) {
256  l = 2 * j ;
257  }
258  else
259  l = 2 * j + 1 ;
260 
261  for (int i=0; i<k_div; i++) {
262  for (int j2=0; j2<k_div; j2++) {
263  matrix.set(i, j2) = cf_cil(l, nr-1-i, j2) ;
264  }
265  }
266 
267  matrix.set_band(k_div, k_div) ;
268 
269  matrix.set_lu() ;
270 
271  for (int i=0; i<k_div; i++) {
272  rhs_cos.set(i) = nilm_cos(m, l, nr-1-i) ;
273  rhs_sin.set(i) = nilm_sin(m, l, nr-1-i) ;
274  }
275 
276  Tbl sol_cos = matrix.inverse(rhs_cos) ;
277  Tbl sol_sin = matrix.inverse(rhs_sin) ;
278 
279  for (int i=0; i<k_div; i++) {
280  alm_cos.set(m, l, i) = sol_cos(i) ;
281  alm_sin.set(m, l, i) = sol_sin(i) ;
282  }
283  }
284  }
285 
286  // -------------------------------------------------------
287  // Construction of the diverging analytic source
288  // -------------------------------------------------------
289 
290  source_div.set_etat_qcq() ;
291  (source_div.va).set_etat_cf_qcq() ;
292  (source_div.va).c_cf->set_etat_qcq() ;
293  (source_div.va).c_cf->t[0]->set_etat_qcq() ;
294 
295  Valeur& sdva = source_div.va ;
296  Mtbl_cf& sdva_cf = *(sdva.c_cf) ;
297 
298  // Initialization
299  for (int k=0; k<=np; k+=2) {
300  for (int j=0; j<nt; j++) {
301  for (int i=0; i<nr; i++) {
302  sdva_cf.set(0, k, j, i) = 0 ;
303  sdva_cf.set(0, k+1, j, i) = 0 ;
304  }
305  }
306  }
307 
308  for (int k=0; k<=np; k+=2) {
309 
310  int m = k / 2 ;
311 
312  for (int j=0; j<nt; j++) {
313 
314  int l ;
315 
316  if (m%2 == 0) {
317  l = 2 * j ;
318  }
319  else
320  l = 2 * j + 1 ;
321 
322  for (int i=0; i<nr; i++) {
323 
324  for (int k_deg=1; k_deg<=k_div; k_deg++) {
325 
326  sdva_cf.set(0, k, j, i) = sdva_cf(0, k, j, i)
327  + alm_cos(m, l, k_deg-1) * cf_cil(l, i, k_deg-1) ;
328  sdva_cf.set(0, k+1, j, i) = sdva_cf(0, k+1, j, i)
329  + alm_sin(m, l, k_deg-1) * cf_cil(l, i, k_deg-1) ;
330 
331  }
332  }
333  }
334  }
335 
336  source_div.annule(nzet, nzm1) ;
337 
338  Base_val base_std = mg->std_base_scal() ;
339 
340  Base_val& base_s_div = sdva.base ;
341  for (int l=0; l<=nzm1; l++) {
342  int base_s_div_r = base_std.b[l] & MSQ_R ;
343  int base_s_div_p = base_std.b[l] & MSQ_P ;
344 
345  base_s_div.b[l] = base_s_div_r | T_LEG_P | base_s_div_p ;
346  }
347 
348  sdva_cf.base = base_s_div ; // copy of the base in the Mtbl_cf
349 
350  sdva.ylm_i() ;
351 
352  // ------------------------------------------------
353  // Construction of the regularized source
354  // ------------------------------------------------
355 
356  source_regu.set_etat_qcq() ;
357  source_regu = source - source_div ;
358 
359  // -------------------------------------------------------------
360  // Solving the Poisson equation for regularized source
361  // -------------------------------------------------------------
362 
363  source_regu.set_dzpuis(4) ;
364 
365  assert(uu_regu.get_mp()->get_mg() == mg) ;
366 
367  (*this).poisson(source_regu, par, uu_regu) ;
368 
369  // ----------------------------------------------------------
370  // Construction of the diverging analytic potential
371  // ----------------------------------------------------------
372 
373  Tbl cf_pil(2*nt, nr, k_div) ;
374  cf_pil.set_etat_qcq() ;
375 
376  double* tmp2 = new double[nr] ;
377 
378  for (int k_deg=1; k_deg<=k_div; k_deg++) {
379 
380  for (int l=0; l<2*nt; l++) {
381 
382  for (int i=0; i<nr; i++) {
383 
384  double xi = grid.x[i] ;
385  tmp2[i] = pow(xi, l) * pow(1. - xi * xi, aa + 1. + k_deg) ;
386 
387  }
388 
389  if (l%2 == 0) {
390  cfrchebp(deg, dim, tmp2, dim, tmp2) ;
391  }
392  else
393  cfrchebi(deg, dim, tmp2, dim, tmp2) ;
394 
395  for (int i=0; i<nr; i++) {
396 
397  cf_pil.set(l, i, k_deg-1) = tmp2[i] ;
398 
399  }
400  }
401  }
402 
403  uu_div.set_etat_qcq() ;
404  (uu_div.va).set_etat_cf_qcq() ;
405  ((uu_div.va).c_cf)->set_etat_qcq() ;
406  ((uu_div.va).c_cf)->t[0]->set_etat_qcq() ;
407 
408  Valeur& udva = uu_div.va ;
409  Mtbl_cf& udva_cf = *(udva.c_cf) ;
410 
411  // Initialization
412  for (int k=0; k<=np; k+=2) {
413  for (int j=0; j<nt; j++) {
414  for (int i=0; i<nr; i++) {
415  udva_cf.set(0, k, j, i) = 0 ;
416  udva_cf.set(0, k+1, j, i) = 0 ;
417  }
418  }
419  }
420 
421  for (int k=0; k<=np; k+=2) {
422 
423  int m = k / 2 ;
424 
425  for (int j=0; j<nt; j++) {
426 
427  int l ;
428 
429  if (m%2 == 0) {
430  l = 2 * j ;
431  }
432  else
433  l = 2 * j + 1 ;
434 
435  for (int i=0; i<nr; i++) {
436 
437  for (int k_deg=1; k_deg<=k_div; k_deg++) {
438 
439  udva_cf.set(0, k, j, i) = udva_cf(0, k, j, i)
440  + alm_cos(m, l, k_deg-1) * cf_pil(l, i, k_deg-1) ;
441  udva_cf.set(0, k+1, j, i) = udva_cf(0, k+1, j, i)
442  + alm_sin(m, l, k_deg-1) * cf_pil(l, i, k_deg-1) ;
443 
444  }
445  }
446  }
447  }
448 
449  uu_div.annule(nzet, nzm1) ;
450 
451  Base_val& base_uu_div = (uu_div.va).base ;
452  for (int l=0; l<=nzm1; l++) {
453  int base_uu_r = base_std.b[l] & MSQ_R ;
454  int base_uu_p = base_std.b[l] & MSQ_P ;
455 
456  base_uu_div.b[l] = base_uu_r | T_LEG_P | base_uu_p ;
457  }
458 
459  udva_cf.base = base_uu_div ; // copy of the base in the Mtbl_cf
460 
461  udva.ylm_i() ;
462 
463  // Changing the radial coordinate from "xi" to "r"
464  // -----------------------------------------------
465 
466  udva = udva * alpha[0] * alpha[0] ;
467 
468  // ---------------------------------------------
469  // Construction of the total potential
470  // ---------------------------------------------
471 
472  uu.set_etat_qcq() ;
473  uu = uu_regu + uu_div ;
474 
475  // -------------------------------------------------------------------
476  // Construction of the derivative of the diverging potential
477  // -------------------------------------------------------------------
478 
479  duu_div.set_etat_qcq() ;
480 
481  duu_div.set(0).set_etat_qcq() ;
482  (duu_div.set(0).va).set_etat_cf_qcq() ;
483  ((duu_div.set(0).va).c_cf)->set_etat_qcq() ;
484  ((duu_div.set(0).va).c_cf)->t[0]->set_etat_qcq() ;
485 
486  duu_div.set(1).set_etat_qcq() ;
487  (duu_div.set(1).va).set_etat_cf_qcq() ;
488  ((duu_div.set(1).va).c_cf)->set_etat_qcq() ;
489  ((duu_div.set(1).va).c_cf)->t[0]->set_etat_qcq() ;
490 
491  duu_div.set(2).set_etat_qcq() ;
492  (duu_div.set(2).va).set_etat_cf_qcq() ;
493  ((duu_div.set(2).va).c_cf)->set_etat_qcq() ;
494  ((duu_div.set(2).va).c_cf)->t[0]->set_etat_qcq() ;
495 
496  Valeur& vr = duu_div.set(0).va ;
497  Valeur& vt = duu_div.set(1).va ;
498  Valeur& vp = duu_div.set(2).va ;
499 
500  Mtbl_cf& vr_cf = *(vr.c_cf) ;
501  Mtbl_cf& vt_cf = *(vt.c_cf) ;
502  Mtbl_cf& vp_cf = *(vp.c_cf) ;
503 
504  // -----------
505  // Radial part
506  // -----------
507 
508  Tbl cf_dril(2*nt, nr, k_div) ;
509  cf_dril.set_etat_qcq() ;
510 
511  double* tmp3 = new double[nr] ;
512 
513  for (int k_deg=1; k_deg<=k_div; k_deg++) {
514 
515  for (int i=0; i<nr; i++) {
516 
517  double xi = grid.x[i] ;
518  tmp3[i] = -2. * (aa + 1. + k_deg) * xi
519  * pow(1. - xi * xi, aa + k_deg) ;
520 
521  }
522 
523  cfrchebi(deg, dim, tmp3, dim, tmp3) ;
524 
525  for (int i=0; i<nr; i++) {
526 
527  cf_dril.set(0, i, k_deg-1) = tmp3[i] ;
528 
529  }
530 
531  for (int l=1; l<2*nt; l++) {
532 
533  for (int i=0; i<nr; i++) {
534 
535  double xi = grid.x[i] ;
536  tmp3[i] = l * pow(xi, l - 1.) * pow(1. - xi * xi, aa + 1. + k_deg)
537  -2. * (aa + 1. + k_deg) * pow(xi, l + 1.)
538  * pow(1. - xi * xi, aa + k_deg) ;
539 
540  }
541 
542  if (l%2 == 0) {
543  cfrchebi(deg, dim, tmp3, dim, tmp3) ;
544  }
545  else
546  cfrchebp(deg, dim, tmp3, dim, tmp3) ;
547 
548  for (int i=0; i<nr; i++) {
549 
550  cf_dril.set(l, i, k_deg-1) = tmp3[i] ;
551 
552  }
553  }
554  }
555 
556  // Initialization
557  for (int k=0; k<=np; k+=2) {
558  for (int j=0; j<nt; j++) {
559  for (int i=0; i<nr; i++) {
560  vr_cf.set(0, k, j, i) = 0 ;
561  vr_cf.set(0, k+1, j, i) = 0 ;
562  }
563  }
564  }
565 
566  for (int k=0; k<=np; k+=2) {
567 
568  int m = k / 2 ;
569 
570  for (int j=0; j<nt; j++) {
571 
572  int l ;
573 
574  if (m%2 == 0) {
575  l = 2 * j ;
576  }
577  else
578  l = 2 * j + 1 ;
579 
580  for (int i=0; i<nr; i++) {
581 
582  for (int k_deg=1; k_deg<=k_div; k_deg++) {
583 
584  vr_cf.set(0, k, j, i) = vr_cf(0, k, j, i)
585  + alm_cos(m, l, k_deg-1) * cf_dril(l, i, k_deg-1) ;
586  vr_cf.set(0, k+1, j, i) = vr_cf(0, k+1, j, i)
587  + alm_sin(m, l, k_deg-1) * cf_dril(l, i, k_deg-1) ;
588 
589  }
590  }
591  }
592  }
593 
594  (duu_div.set(0)).annule(nzet, nzm1) ;
595 
596  // Reconstruction of the basis of the radial part
597  // ----------------------------------------------
598 
599  Base_val& base_duu_div_r = vr.base ;
600  for (int l=0; l<=nzm1; l++) {
601  int base_duu_r_p = base_std.b[l] & MSQ_P ;
602 
603  base_duu_div_r.b[l] = R_CHEBPIM_I | T_LEG_P | base_duu_r_p ;
604  }
605 
606  vr_cf.base = base_duu_div_r ;
607  vr.ylm_i() ;
608 
609  const Coord& RR = dxdr ;
610 
611  // Changing the radial coordinate from "xi" to "r"
612  // -----------------------------------------------
613 
614  Base_val sauve_base( vr.base ) ;
615  vr = duu_div(0).va * alpha[0] * alpha[0] * RR ;
616  vr.base = sauve_base ;
617 
618  // -------------------------
619  // Polar and azimuthal parts
620  // -------------------------
621 
622  Tbl cf_dpil(2*nt, nr, k_div) ;
623  cf_dpil.set_etat_qcq() ;
624 
625  double* tmp4 = new double[nr] ;
626 
627  for (int k_deg=1; k_deg<=k_div; k_deg++) {
628 
629  for (int i=0; i<nr; i++) {
630  tmp4[i] = 0 ;
631  }
632 
633  cfrchebi(deg, dim, tmp4, dim, tmp4) ;
634 
635  for (int i=0; i<nr; i++) {
636 
637  cf_dpil.set(0, i, k_deg-1) = tmp4[i] ;
638 
639  }
640 
641  for (int l=1; l<2*nt; l++) {
642 
643  for (int i=0; i<nr; i++) {
644 
645  double xi = grid.x[i] ;
646  tmp4[i] = pow(xi, l - 1.) * pow(1. - xi * xi, aa + 1. + k_deg) ;
647 
648  }
649 
650  if (l%2 == 0) {
651  cfrchebi(deg, dim, tmp4, dim, tmp4) ;
652  }
653  else
654  cfrchebp(deg, dim, tmp4, dim, tmp4) ;
655 
656  for (int i=0; i<nr; i++) {
657 
658  cf_dpil.set(l, i, k_deg-1) = tmp4[i] ;
659 
660  }
661  }
662  }
663 
664  // Initialization
665  for (int k=0; k<=np; k+=2) {
666  for (int j=0; j<nt; j++) {
667  for (int i=0; i<nr; i++) {
668  vt_cf.set(0, k, j, i) = 0 ;
669  vt_cf.set(0, k+1, j, i) = 0 ;
670  vp_cf.set(0, k, j, i) = 0 ;
671  vp_cf.set(0, k+1, j, i) = 0 ;
672  }
673  }
674  }
675 
676  for (int k=0; k<=np; k+=2) {
677 
678  int m = k / 2 ;
679 
680  for (int j=0; j<nt; j++) {
681 
682  int l ;
683 
684  if (m%2 == 0) {
685  l = 2 * j ;
686  }
687  else
688  l = 2 * j + 1 ;
689 
690  for (int i=0; i<nr; i++) {
691 
692  for (int k_deg=1; k_deg<=k_div; k_deg++) {
693 
694  vt_cf.set(0, k, j, i) = vt_cf(0, k, j, i)
695  + alm_cos(m, l, k_deg-1) * cf_dpil(l, i, k_deg-1) ;
696  vt_cf.set(0, k+1, j, i) = vt_cf(0, k+1, j, i)
697  + alm_sin(m, l, k_deg-1) * cf_dpil(l, i, k_deg-1) ;
698 
699  vp_cf.set(0, k, j, i) = vp_cf(0, k, j, i)
700  + alm_cos(m, l, k_deg-1) * cf_dpil(l, i, k_deg-1) ;
701  vp_cf.set(0, k+1, j, i) = vp_cf(0, k+1, j, i)
702  + alm_sin(m, l, k_deg-1) * cf_dpil(l, i, k_deg-1) ;
703 
704  }
705  }
706  }
707  }
708 
709  (duu_div.set(1)).annule(nzet, nzm1) ;
710  (duu_div.set(2)).annule(nzet, nzm1) ;
711 
712 
713  // Reconstruction of the basis of the polar part
714  // ---------------------------------------------
715 
716  Base_val& base_duu_div_p = vt.base ;
717  for (int l=0; l<=nzm1; l++) {
718  int base_duu_p_p = base_std.b[l] & MSQ_P ;
719 
720  base_duu_div_p.b[l] = R_CHEBPIM_I | T_LEG_P | base_duu_p_p ;
721  }
722 
723  vt_cf.base = base_duu_div_p ;
724  vt.ylm_i() ;
725 
726 
727  // Reconstruction of the basis of the azimuthal part
728  // -------------------------------------------------
729 
730  Base_val& base_duu_div_t = vp.base ;
731  for (int l=0; l<=nzm1; l++) {
732  int base_duu_t_p = base_std.b[l] & MSQ_P ;
733 
734  base_duu_div_t.b[l] = R_CHEBPIM_I | T_LEG_P | base_duu_t_p ;
735  }
736 
737  vp_cf.base = base_duu_div_t ;
738  vp.ylm_i() ;
739 
740 
741  // Calculation of the derivatives
742  // ------------------------------
743 
744  vt = (duu_div(1).va).dsdt() ;
745 
746  vp = (duu_div(2).va).stdsdp() ;
747 
748  // Changing the radial coordinate from "xi" to "r"
749  // -----------------------------------------------
750 
751  vt = duu_div(1).va * alpha[0] ;
752 
753  vp = duu_div(2).va * alpha[0] ;
754 
755  // Set the basis of duu_div to the spherical one
756  // ---------------------------------------------
757 
758  duu_div.set_triad( (*this).get_bvect_spher() ) ;
759 
760  delete [] tmp1 ;
761  delete [] tmp2 ;
762  delete [] tmp3 ;
763  delete [] tmp4 ;
764 
765 
766 }
767 }
Bases of the spectral expansions.
Definition: base_val.h:322
int * b
Array (size: nzone ) of the spectral basis in each domain.
Definition: base_val.h:331
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
Valeur va
The numerical value of the Cmp
Definition: cmp.h:464
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:348
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:304
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
void set_dzpuis(int)
Set a value to dzpuis.
Definition: cmp.C:654
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
3D grid class in one domain.
Definition: grilles.h:194
double * x
Array of values of at the nr collocation points.
Definition: grilles.h:209
virtual void poisson_regular(const Cmp &source, int k_div, int nzet, double unsgam1, Param &par, Cmp &uu, Cmp &uu_regu, Cmp &uu_div, Tenseur &duu_div, Cmp &source_regu, Cmp &source_div) const
Computes the solution of a scalar Poisson equation.
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition: map.h:2033
virtual void dsdt(const Scalar &uu, Scalar &resu) const
Computes of a Scalar.
Definition: map_af_deriv.C:797
virtual void stdsdp(const Scalar &uu, Scalar &resu) const
Computes of a Scalar.
Definition: map_af_deriv.C:823
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1560
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition: map.h:676
Matrix handling.
Definition: matrice.h:152
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:175
double & set(int j, int i)
Read/write of a particuliar element.
Definition: matrice.h:277
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition: matrice.C:424
void set_band(int up, int low) const
Calculate the band storage of *std.
Definition: matrice.C:364
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition: matrice.C:392
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
Definition: grilles.h:500
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Base_val std_base_scal() const
Returns the standard spectral bases for a scalar.
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:186
Base_val base
Bases of the spectral expansions.
Definition: mtbl_cf.h:200
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:294
Parameter storage.
Definition: param.h:125
Basic array class.
Definition: tbl.h:161
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition: tenseur.C:674
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
int get_etat() const
Returns the logical state.
Definition: valeur.h:726
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
Definition: valeur.h:729
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
#define MSQ_R
Extraction de l'info sur R.
Definition: type_parite.h:152
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
Definition: type_parite.h:178
#define T_LEG_P
fct. de Legendre associees paires
Definition: type_parite.h:216
#define MSQ_P
Extraction de l'info sur Phi.
Definition: type_parite.h:156
Lorene prototypes.
Definition: app_hor.h:64