LORENE
scalar_asymptot.C
1/*
2 * Function Scalar::asymptot
3 *
4 */
5
6/*
7 * Copyright (c) 2003 Eric Gourgoulhon & Jerome Novak
8 *
9 * Copyright (c) 1999-2002 Eric Gourgoulhon (for previous class Cmp)
10 * Copyright (c) 1999-2001 Philippe Grandclement (for previous class Cmp)
11 *
12 * This file is part of LORENE.
13 *
14 * LORENE is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU General Public License as published by
16 * the Free Software Foundation; either version 2 of the License, or
17 * (at your option) any later version.
18 *
19 * LORENE is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with LORENE; if not, write to the Free Software
26 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27 *
28 */
29
30
31char scalar_asymptot_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_asymptot.C,v 1.5 2014/10/13 08:53:46 j_novak Exp $" ;
32
33/*
34 * $Id: scalar_asymptot.C,v 1.5 2014/10/13 08:53:46 j_novak Exp $
35 * $Log: scalar_asymptot.C,v $
36 * Revision 1.5 2014/10/13 08:53:46 j_novak
37 * Lorene classes and functions now belong to the namespace Lorene.
38 *
39 * Revision 1.4 2014/10/06 15:16:15 j_novak
40 * Modified #include directives to use c++ syntax.
41 *
42 * Revision 1.3 2004/02/21 17:05:14 e_gourgoulhon
43 * Method Scalar::point renamed Scalar::val_grid_point.
44 * Method Scalar::set_point renamed Scalar::set_grid_point.
45 *
46 * Revision 1.2 2003/10/08 14:24:09 j_novak
47 * replaced mult_r_zec with mult_r_ced
48 *
49 * Revision 1.1 2003/09/25 07:18:00 j_novak
50 * Method asymptot implemented.
51 *
52 * $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_asymptot.C,v 1.5 2014/10/13 08:53:46 j_novak Exp $
53 *
54 */
55
56// Headers C
57#include <cmath>
58
59// Headers Lorene
60#include "tensor.h"
61
62namespace Lorene {
63Valeur** Scalar::asymptot(int n0, const int flag) const {
64
65 assert(n0 >= 0) ;
66 const Mg3d& mg = *(mp->get_mg()) ;
67 int nz = mg.get_nzone() ;
68 int nzm1 = nz-1 ;
69 assert(mg.get_type_r(nzm1) == UNSURR) ;
70 int np = mg.get_np(nzm1) ;
71 int nt = mg.get_nt(nzm1) ;
72 int nr = mg.get_nr(nzm1) ;
73 int nrm1 = nr-1 ;
74
75 Valeur** vu = new Valeur*[n0+1] ;
76 for (int h=0; h<=n0; h++) {
77 vu[h] = new Valeur(mg.get_angu()) ;
78 }
79
80 Scalar uu = *this ;
81
82 int precis = 2 ;
83
84 // The terms 1/r^h with h < dzpuis are null :
85 // -----------------------------------------
86 for (int h=0; h<dzpuis; h++) {
87
88 vu[h]->set_etat_zero() ;
89
90 }
91
92 // Terms 1/r^h with h >= dzpuis :
93 // -----------------------------
94 for (int h=dzpuis; h<=n0; h++) {
95
96 // Value on the sphere S^2 at infinity
97 // -----------------------------------
98 vu[h]->set_etat_c_qcq() ;
99 vu[h]->c->set_etat_qcq() ;
100 for (int l=0; l<nzm1; l++) {
101 vu[h]->c->t[l]->set_etat_zero() ;
102 }
103 vu[h]->c->t[nzm1]->set_etat_qcq() ;
104
105 for (int k=0; k<np; k++) {
106 for (int j=0; j<nt; j++) {
107 vu[h]->set(nzm1, k, j, 0) = uu.val_grid_point(nzm1, k, j, nrm1) ;
108 }
109 }
110
111 vu[h]->set_base( uu.va.base ) ;
112
113 // Printing
114 // --------
115 if (flag != 0) {
116 cout << "Term in 1/r^" << h << endl ;
117 cout << "-------------" << endl ;
118
119 double ndec = 0 ;
120 double vmin = (*vu[h])(nzm1, 0, 0, 0) ;
121 double vmax = vmin ;
122
123 cout << " Values at the point (phi_k, theta_j) : " << endl ;
124 cout.precision(precis) ;
125 cout.setf(ios::showpoint) ;
126 for (int k=0; k<np; k++) {
127 cout << " k=" << k << " : " ;
128 for (int j=0; j<nt; j++) {
129 double xx = (*vu[h])(nzm1, k, j, 0) ;
130 cout << " " << setw(precis) << xx ;
131 ndec += fabs(xx) ;
132 vmin = ( xx < vmin ) ? xx : vmin ;
133 vmax = ( xx > vmax ) ? xx : vmax ;
134 }
135 cout << endl;
136 }
137 ndec /= np*nt ;
138 cout << "Minimum value on S^2 : " << vmin << endl ;
139 cout << "Maximum value on S^2 : " << vmax << endl ;
140 cout << "L^1 norm on S^2 : " << ndec << endl ;
141 }
142 // The value at infinity is substracted
143 // ------------------------------------
144 for (int k=0; k<np; k++) {
145 for (int j=0; j<nt; j++) {
146 double v_inf = (*vu[h])(nzm1, k, j, 0) ;
147 for (int i=0; i<nr; i++) {
148 uu.set_grid_point(nzm1, k, j, i) -= v_inf ;
149 }
150 }
151 }
152
153 // Mutliplication by r
154 // -------------------
155
156 uu.mult_r_ced() ;
157
158 } // End of loop on h (development order)
159
160 return vu ;
161
162}
163}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Multi-domain grid.
Definition grilles.h:273
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition mg3d.C:473
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
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:474
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Valeur ** asymptot(int n, const int flag=0) const
Asymptotic expansion at r = infinity.
int dzpuis
Power of r by which the quantity represented by this must be divided in the compactified external d...
Definition scalar.h:403
Values and coefficients of a (real-value) function.
Definition valeur.h:287
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Lorene prototypes.
Definition app_hor.h:64