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
28char 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
107namespace Lorene {
108
109void 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
void set_dzpuis(int)
Set a value to dzpuis.
Definition cmp.C:654
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
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.
virtual void stdsdp(const Scalar &uu, Scalar &resu) const
Computes of a Scalar.
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
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_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
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
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
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
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
void coef() const
Computes the coeffcients of *this.
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
Definition valeur.h:729
void ylm_i()
Inverse of ylm()
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.
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
#define T_LEG_P
fct. de Legendre associees paires
#define MSQ_P
Extraction de l'info sur Phi.
Lorene prototypes.
Definition app_hor.h:64