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