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 
28 char 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
74 namespace Lorene {
75 void 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
180 void 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 
283 void Map_et::lapang(const Scalar& , Scalar& ) const {
284 
285  cout << "Map_et::lapang : not implemented yet!" << endl ;
286  abort() ;
287 
288 }
289 
290 
291 }
Lorene::Base_val
Bases of the spectral expansions.
Definition: base_val.h:322
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
Lorene::Cmp::check_dzpuis
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
Lorene::Map_radial::d2rdtdx
Coord d2rdtdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1640
Lorene::Valeur::stdsdp
const Valeur & stdsdp() const
Returns of *this.
Definition: valeur_stdsdp.C:60
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Map_radial::xsr
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1549
Lorene::Mtbl
Multi-domain array.
Definition: mtbl.h:118
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp
Definition: cmp.h:464
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Cmp::get_etat
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
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:474
Lorene::Scalar::check_dzpuis
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
Lorene::Valeur::annule
void annule(int l)
Sets the Valeur to zero in a given domain.
Definition: valeur.C:744
Lorene::Cmp::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: cmp.C:289
Lorene::Map_radial::sstd2rdpdx
Coord sstd2rdpdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1648
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:861
Lorene::Map_et::laplacien
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
Lorene::Scalar::get_spectral_va
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:601
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
Lorene::Map_radial::sr2drdt
Coord sr2drdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1600
Lorene::Map_radial::dxdr
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1560
Lorene::Map_radial::sr2stdrdp
Coord sr2stdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1608
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Lorene::Map_radial::srstdrdp
Coord srstdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1592
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Valeur::set_base
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition: valeur.C:810
Lorene::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:353
Lorene::Map_et::lapang
virtual void lapang(const Scalar &uu, Scalar &lap) const
Computes the angular Laplacian of a scalar field.
Definition: map_et_lap.C:283
Lorene::Scalar::set_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:324
Lorene::Scalar::set_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:808
Lorene::Map_radial::d2rdx2
Coord d2rdx2
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1619
Lorene::Valeur::d2sdx2
const Valeur & d2sdx2() const
Returns of *this.
Definition: valeur_d2sdx2.C:114
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
Lorene::Valeur::dsdt
const Valeur & dsdt() const
Returns of *this.
Definition: valeur_dsdt.C:112
Lorene::Valeur::dsdx
const Valeur & dsdx() const
Returns of *this.
Definition: valeur_dsdx.C:111
Lorene::Valeur::lapang
const Valeur & lapang() const
Returns the angular Laplacian of *this.
Definition: valeur_lapang.C:72
Lorene::Valeur::base
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
Lorene::Valeur::sx
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR )
Definition: valeur_sx.C:110
Lorene::Map_radial::srdrdt
Coord srdrdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1584
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
Lorene::Cmp::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:304
Lorene::Map_radial::lapr_tp
Coord lapr_tp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1631
Lorene::Cmp::get_mp
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
Lorene::Cmp::set_dzpuis
void set_dzpuis(int)
Set a value to dzpuis.
Definition: cmp.C:654
Lorene::Map::mg
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition: map.h:676