LORENE
map_et_lap.C
1/*
2 * Computes the Laplacian of a scalar field (represented by a Cmp) when
3 * the mapping belongs to the Map_et class
4 */
5
6/*
7 * Copyright (c) 1999-2001 Eric Gourgoulhon
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_et_lap_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_et_lap.C,v 1.4 2014/10/13 08:53:05 j_novak Exp $" ;
29
30/*
31 * $Id: map_et_lap.C,v 1.4 2014/10/13 08:53:05 j_novak Exp $
32 * $Log: map_et_lap.C,v $
33 * Revision 1.4 2014/10/13 08:53:05 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.3 2005/11/24 09:25:07 j_novak
37 * Added the Scalar version for the Laplacian
38 *
39 * Revision 1.2 2003/10/15 16:03:37 j_novak
40 * Added the angular Laplace operator for Scalar.
41 *
42 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
43 * LORENE
44 *
45 * Revision 1.4 2000/02/25 09:02:18 eric
46 * Remplacement de (uu.get_dzpuis() == 0) par uu.check_dzpuis(0).
47 *
48 * Revision 1.3 2000/01/26 13:11:07 eric
49 * Reprototypage complet :
50 * le resultat est desormais suppose alloue a l'exterieur de la routine
51 * et est passe en argument (Cmp& resu).
52 * .
53 *
54 * Revision 1.2 2000/01/14 14:55:05 eric
55 * Suppression de l'assert(sauve_base == vresu.base)
56 * car sauve_base == vresu.base n'est pas necessairement vrai (cela
57 * depend de l'histoire du Cmp uu, notamment de si uu.p_dsdx est
58 * a jour).
59 *
60 * Revision 1.1 1999/12/20 17:27:30 eric
61 * Initial revision
62 *
63 *
64 * $Header: /cvsroot/Lorene/C++/Source/Map/map_et_lap.C,v 1.4 2014/10/13 08:53:05 j_novak Exp $
65 *
66 */
67
68// Header Lorene
69#include "cmp.h"
70#include "tensor.h"
71
72
73// Laplacian: Scalar version
74namespace Lorene {
75void Map_et::laplacien(const Scalar& uu, int zec_mult_r, Scalar& resu) const {
76
77 assert (uu.get_etat() != ETATNONDEF) ;
78 assert (uu.get_mp().get_mg() == mg) ;
79
80 // Particular case of null value:
81
82 if ((uu.get_etat() == ETATZERO) || (uu.get_etat() == ETATUN)) {
83 resu.set_etat_zero() ;
84 return ;
85 }
86
87 assert( uu.get_etat() == ETATQCQ ) ;
88 assert( uu.check_dzpuis(0) ) ;
89
90 int nz = get_mg()->get_nzone() ;
91 int nzm1 = nz - 1 ;
92
93 // Indicator of existence of a compactified external domain
94 bool zec = false ;
95 if (mg->get_type_r(nzm1) == UNSURR) {
96 zec = true ;
97 }
98
99 if ( zec && (zec_mult_r != 4) ) {
100 cout << "Map_et::laplacien : the case zec_mult_r = " <<
101 zec_mult_r << " is not implemented !" << endl ;
102 abort() ;
103 }
104
105 //--------------------
106 // First operations
107 //--------------------
108
109 Valeur duudx = uu.get_spectral_va().dsdx() ; // d/dx
110 Valeur d2uudx2 = uu.get_spectral_va().d2sdx2() ; // d^2/dx^2
111
112 const Valeur& d2uudtdx = duudx.dsdt() ; // d^2/dxdtheta
113
114 const Valeur& std2uudpdx = duudx.stdsdp() ; // 1/sin(theta) d^2/dxdphi
115
116 //------------------
117 // Angular Laplacian
118 //------------------
119
120 Valeur sxlapang = uu.get_spectral_va() ;
121
122 sxlapang.ylm() ;
123
124 sxlapang = sxlapang.lapang() ;
125
126 sxlapang = sxlapang.sx() ; // 1/x in the nucleus
127 // Id in the shells
128 // 1/(x-1) in the ZEC
129
130 //------------------------------------
131 // (2 dx/dR d/dx + x/R 1/x Lap_ang)/x
132 //------------------------------------
133
134 Valeur varduudx = duudx ;
135
136 if (zec) {
137 varduudx.annule(nzm1) ; // term in d/dx set to zero in the ZEC
138 }
139
140 resu.set_etat_qcq() ;
141
142 Valeur& vresu = resu.set_spectral_va() ;
143
144 Base_val sauve_base = varduudx.base ;
145
146 vresu = double(2) * dxdr * varduudx + xsr * sxlapang ;
147
148 vresu.set_base(sauve_base) ;
149
150 vresu = vresu.sx() ;
151
152 //--------------
153 // Final result
154 //--------------
155
156 Mtbl unjj = double(1) + srdrdt*srdrdt + srstdrdp*srstdrdp ;
157
158 sauve_base = d2uudx2.base ;
159// assert(sauve_base == vresu.base) ; // this is not necessary true
160
161 vresu = dxdr*dxdr * unjj * d2uudx2 + xsr * vresu
162 - double(2) * dxdr * ( sr2drdt * d2uudtdx
163 + sr2stdrdp * std2uudpdx ) ;
164
165 vresu += - dxdr * ( lapr_tp + dxdr * (
166 dxdr* unjj * d2rdx2
167 - double(2) * ( sr2drdt * d2rdtdx + sr2stdrdp * sstd2rdpdx ) )
168 ) * duudx ;
169
170 vresu.set_base(sauve_base) ;
171
172 if (zec == 1) {
173 resu.set_dzpuis(zec_mult_r) ;
174 }
175
176}
177
178
179// Laplacian: Cmp version
180void Map_et::laplacien(const Cmp& uu, int zec_mult_r, Cmp& resu) const {
181
182 assert (uu.get_etat() != ETATNONDEF) ;
183 assert (uu.get_mp()->get_mg() == mg) ;
184
185 // Particular case of null value:
186
187 if (uu.get_etat() == ETATZERO) {
188 resu.set_etat_zero() ;
189 return ;
190 }
191
192 assert( uu.get_etat() == ETATQCQ ) ;
193 assert( uu.check_dzpuis(0) ) ;
194
195 int nz = get_mg()->get_nzone() ;
196 int nzm1 = nz - 1 ;
197
198 // Indicator of existence of a compactified external domain
199 bool zec = false ;
200 if (mg->get_type_r(nzm1) == UNSURR) {
201 zec = true ;
202 }
203
204 if ( zec && (zec_mult_r != 4) ) {
205 cout << "Map_et::laplacien : the case zec_mult_r = " <<
206 zec_mult_r << " is not implemented !" << endl ;
207 abort() ;
208 }
209
210 //--------------------
211 // First operations
212 //--------------------
213
214 Valeur duudx = (uu.va).dsdx() ; // d/dx
215 Valeur d2uudx2 = (uu.va).d2sdx2() ; // d^2/dx^2
216
217 const Valeur& d2uudtdx = duudx.dsdt() ; // d^2/dxdtheta
218
219 const Valeur& std2uudpdx = duudx.stdsdp() ; // 1/sin(theta) d^2/dxdphi
220
221 //------------------
222 // Angular Laplacian
223 //------------------
224
225 Valeur sxlapang = uu.va ;
226
227 sxlapang.ylm() ;
228
229 sxlapang = sxlapang.lapang() ;
230
231 sxlapang = sxlapang.sx() ; // 1/x in the nucleus
232 // Id in the shells
233 // 1/(x-1) in the ZEC
234
235 //------------------------------------
236 // (2 dx/dR d/dx + x/R 1/x Lap_ang)/x
237 //------------------------------------
238
239 Valeur varduudx = duudx ;
240
241 if (zec) {
242 varduudx.annule(nzm1) ; // term in d/dx set to zero in the ZEC
243 }
244
245 resu.set_etat_qcq() ;
246
247 Valeur& vresu = resu.va ;
248
249 Base_val sauve_base = varduudx.base ;
250
251 vresu = double(2) * dxdr * varduudx + xsr * sxlapang ;
252
253 vresu.set_base(sauve_base) ;
254
255 vresu = vresu.sx() ;
256
257 //--------------
258 // Final result
259 //--------------
260
261 Mtbl unjj = double(1) + srdrdt*srdrdt + srstdrdp*srstdrdp ;
262
263 sauve_base = d2uudx2.base ;
264// assert(sauve_base == vresu.base) ; // this is not necessary true
265
266 vresu = dxdr*dxdr * unjj * d2uudx2 + xsr * vresu
267 - double(2) * dxdr * ( sr2drdt * d2uudtdx
268 + sr2stdrdp * std2uudpdx ) ;
269
270 vresu += - dxdr * ( lapr_tp + dxdr * (
271 dxdr* unjj * d2rdx2
272 - double(2) * ( sr2drdt * d2rdtdx + sr2stdrdp * sstd2rdpdx ) )
273 ) * duudx ;
274
275 vresu.set_base(sauve_base) ;
276
277 if (zec == 1) {
278 resu.set_dzpuis(zec_mult_r) ;
279 }
280
281}
282
283void Map_et::lapang(const Scalar& , Scalar& ) const {
284
285 cout << "Map_et::lapang : not implemented yet!" << endl ;
286 abort() ;
287
288}
289
290
291}
Bases of the spectral expansions.
Definition base_val.h:322
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 set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition cmp.C:304
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition cmp.C:289
void set_dzpuis(int)
Set a value to dzpuis.
Definition cmp.C:654
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition cmp.C:715
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
virtual void lapang(const Scalar &uu, Scalar &lap) const
Computes the angular Laplacian of a scalar field.
Definition map_et_lap.C:283
virtual void laplacien(const Scalar &uu, int zec_mult_r, Scalar &lap) const
Computes the Laplacian of a scalar field.
Definition map_et_lap.C:75
Coord d2rdx2
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1619
Coord sr2drdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1600
Coord srstdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1592
Coord d2rdtdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1640
Coord sstd2rdpdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1648
Coord lapr_tp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1631
Coord sr2stdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1608
Coord srdrdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1584
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1549
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
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:474
Multi-domain array.
Definition mtbl.h:118
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition scalar.C:873
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition scalar.h:601
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
Values and coefficients of a (real-value) function.
Definition valeur.h:287
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR )
Definition valeur_sx.C:110
const Valeur & stdsdp() const
Returns of *this.
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition valeur.C:810
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
const Valeur & dsdt() const
Returns of *this.
void annule(int l)
Sets the Valeur to zero in a given domain.
Definition valeur.C:744
const Valeur & dsdx() const
Returns of *this.
const Valeur & d2sdx2() const
Returns of *this.
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
const Valeur & lapang() const
Returns the angular Laplacian of *this.
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene prototypes.
Definition app_hor.h:64