LORENE
scalar_deriv.C
1 /*
2  * Computations of partial derivatives d/dx, d/dy and d/dz of a Scalar.
3  */
4 
5 /*
6  * Copyright (c) 2003-2004 Eric Gourgoulhon & Jerome Novak
7  *
8  * Copyright (c) 1999-2001 Eric Gourgoulhon (for a preceding Cmp version)
9  * Copyright (c) 1999-2001 Philippe Grandclement (for a preceding Cmp version)
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 
32 
33 
34 /*
35  * Revision 1.17 2005/11/24 09:25:08 j_novak
36  * Added the Scalar version for the Laplacian
37  *
38  * Revision 1.16 2005/09/15 15:51:27 j_novak
39  * The "rotation" (change of triad) methods take now Scalars as default
40  * arguments.
41  *
42  * Revision 1.15 2005/05/25 16:11:05 j_novak
43  * Better handling of the case with no compactified domain.
44  *
45  * Revision 1.14 2004/08/24 09:14:52 p_grandclement
46  * Addition of some new operators, like Poisson in 2d... It now requieres the
47  * GSL library to work.
48  *
49  * Also, the way a variable change is stored by a Param_elliptic is changed and
50  * no longer uses Change_var but rather 2 Scalars. The codes using that feature
51  * will requiere some modification. (It should concern only the ones about monopoles)
52  *
53  * Revision 1.13 2004/06/22 08:50:00 p_grandclement
54  * Addition of everything needed for using the logarithmic mapping
55  *
56  * Revision 1.12 2004/02/26 22:51:34 e_gourgoulhon
57  * Added methods derive_cov, derive_con and derive_lie.
58  *
59  * Revision 1.11 2004/01/28 14:02:06 j_novak
60  * Suppressed base handling.
61  *
62  * Revision 1.10 2004/01/28 13:25:42 j_novak
63  * The ced_mult_r arguments have been suppressed from the Scalar::*dsd* methods.
64  * In the div/mult _r_dzpuis, there is no more default value.
65  *
66  * Revision 1.9 2004/01/27 15:10:02 j_novak
67  * New methods Scalar::div_r_dzpuis(int) and Scalar_mult_r_dzpuis(int)
68  * which replace div_r_inc*. Tried to clean the dzpuis handling.
69  * WARNING: no testing at this point!!
70  *
71  * Revision 1.8 2003/11/03 13:37:59 j_novak
72  * Still dzpuis...
73  *
74  * Revision 1.7 2003/10/29 13:14:03 e_gourgoulhon
75  * Added integer argument to derivative functions dsdr, etc...
76  * so that one can choose the dzpuis of the result (default=2).
77  *
78  * Revision 1.6 2003/10/17 13:46:15 j_novak
79  * The argument is now between 1 and 3 (instead of 0->2)
80  *
81  * Revision 1.5 2003/10/15 16:03:38 j_novak
82  * Added the angular Laplace operator for Scalar.
83  *
84  * Revision 1.4 2003/10/15 10:43:58 e_gourgoulhon
85  * Added new methods dsdt and stdsdp.
86  *
87  * Revision 1.3 2003/10/11 14:43:29 e_gourgoulhon
88  * Changed name of local Cmp "deriv" to "derivee" (in order not
89  * to shadow the member deriv).
90  *
91  * Revision 1.2 2003/10/10 15:57:29 j_novak
92  * Added the state one (ETATUN) to the class Scalar
93  *
94  * Revision 1.1 2003/09/25 08:13:52 j_novak
95  * Added method for calculating derivatives
96  *
97  *
98  * $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_deriv.C,v 1.21 2016/12/05 16:18:18 j_novak Exp $
99  *
100  */
101 
102 // Headers Lorene
103 #include "scalar.h"
104 #include "map.h"
105 #include "tensor.h"
106 #include "cmp.h"
107 
108  //---------------------//
109  // d/dr //
110  //---------------------//
111 
112 namespace Lorene {
113 const Scalar& Scalar::dsdr() const {
114 
115  // Protection
116  assert(etat != ETATNONDEF) ;
117 
118  // If the derivative has not been previously computed, the
119  // computation must be done by the appropriate routine of the mapping :
120 
121  if (p_dsdr == 0x0) {
122  p_dsdr = new Scalar(*mp) ;
123  if (etat == ETATUN) {
124  p_dsdr->set_etat_zero() ;
125  }
126  else {
127  mp->dsdr(*this, *p_dsdr) ;
128 
129  }
130  }
131 
132  int dzp = (dzpuis == 0) ? 2 : dzpuis+1 ;
133  if (mp->get_mg()->get_type_r(mp->get_mg()->get_nzone() - 1) != UNSURR)
134  dzp = 0 ;
135  p_dsdr->set_dzpuis(dzp) ;
136 
137  return *p_dsdr ;
138 
139 }
140 
141  //--------------------//
142  // 1/r d/dtheta //
143  //--------------------//
144 
145 const Scalar& Scalar::srdsdt() const {
146 
147  // Protection
148  assert(etat != ETATNONDEF) ;
149 
150  // If the derivative has not been previously computed, the
151  // computation must be done by the appropriate routine of the mapping :
152 
153  if (p_srdsdt == 0x0) {
154  p_srdsdt = new Scalar(*mp) ;
155  if (etat == ETATUN) {
157  }
158  else {
159  mp->srdsdt(*this, *p_srdsdt) ;
160  }
161  }
162 
163  int dzp = (dzpuis == 0) ? 2 : dzpuis+1 ;
164  if (mp->get_mg()->get_type_r(mp->get_mg()->get_nzone() - 1) != UNSURR)
165  dzp = 0 ;
166  p_srdsdt->set_dzpuis(dzp) ;
167 
168  return *p_srdsdt ;
169 
170 }
171 
172 
173  //------------------------------//
174  // 1/(r sin(theta)) d/dphi //
175  //------------------------------//
176 
177 const Scalar& Scalar::srstdsdp() const {
178 
179  // Protection
180  assert(etat != ETATNONDEF) ;
181 
182  // If the derivative has not been previously computed, the
183  // computation must be done by the appropriate routine of the mapping :
184 
185  if (p_srstdsdp == 0x0) {
186  p_srstdsdp = new Scalar(*mp) ;
187  if (etat == ETATUN) {
189  }
190  else {
191  mp->srstdsdp(*this, *p_srstdsdp) ;
192  }
193  }
194 
195  int dzp = (dzpuis == 0) ? 2 : dzpuis+1 ;
196  if (mp->get_mg()->get_type_r(mp->get_mg()->get_nzone() - 1) != UNSURR)
197  dzp = 0 ;
198  p_srstdsdp->set_dzpuis(dzp) ;
199 
200  return *p_srstdsdp ;
201 
202 }
203 
204  //--------------------//
205  // d/dtheta //
206  //--------------------//
207 
208 const Scalar& Scalar::dsdt() const {
209 
210  assert(etat != ETATNONDEF) ; // Protection
211 
212  // If the derivative has not been previously computed, the
213  // computation must be done by the appropriate routine of the mapping :
214 
215  if (p_dsdt == 0x0) {
216 
217  p_dsdt = new Scalar(*mp) ;
218 
219  if (etat == ETATUN) {
220  p_dsdt->set_etat_zero() ;
221  }
222  else {
223  mp->dsdt(*this, *p_dsdt) ;
224  }
225  }
226 
227 
229 
230  return *p_dsdt ;
231 
232 }
233 
234  //------------------------------//
235  // 1/sin(theta) d/dphi //
236  //------------------------------//
237 
238 const Scalar& Scalar::stdsdp() const {
239 
240  assert(etat != ETATNONDEF) ; // Protection
241 
242  // If the derivative has not been previously computed, the
243  // computation must be done by the appropriate routine of the mapping :
244 
245  if (p_stdsdp == 0x0) {
246 
247  p_stdsdp = new Scalar(*mp) ;
248 
249  if (etat == ETATUN) {
250  p_stdsdp->set_etat_zero() ;
251  }
252  else {
253  mp->stdsdp(*this, *p_stdsdp) ;
254  }
255  }
257 
258  return *p_stdsdp ;
259 
260 }
261 
262  //-----------------//
263  // d/dx //
264  //-----------------//
265 
266 const Scalar& Scalar::dsdx() const {
267 
268  // Protection
269  assert(etat != ETATNONDEF) ;
270 
271  // If the derivative has not been previously computed, the
272  // computation must be done by the appropriate routine of the mapping :
273 
274  if (p_dsdx == 0x0) {
275  p_dsdx = new Scalar(*mp) ;
276  if (etat == ETATUN) {
277  p_dsdx->set_etat_zero() ;
278  }
279  else {
281  }
282  }
283 
284  int dzp = (dzpuis == 0) ? 2 : dzpuis+1 ;
285  if (mp->get_mg()->get_type_r(mp->get_mg()->get_nzone() - 1) != UNSURR)
286  dzp = 0 ;
287  p_dsdx->set_dzpuis(dzp) ;
288 
289  return *p_dsdx ;
290 
291 }
292 
293  //-----------------//
294  // d/dy //
295  //-----------------//
296 
297 const Scalar& Scalar::dsdy() const {
298 
299  // Protection
300  assert(etat != ETATNONDEF) ;
301 
302  // If the derivative has not been previously computed, the
303  // computation must be done by the appropriate routine of the mapping :
304 
305  if (p_dsdy == 0x0) {
306  p_dsdy = new Scalar(*mp) ;
307  if (etat == ETATUN) {
308  p_dsdy->set_etat_zero() ;
309  }
310  else {
312  }
313  }
314 
315  int dzp = (dzpuis == 0) ? 2 : dzpuis+1 ;
316  if (mp->get_mg()->get_type_r(mp->get_mg()->get_nzone() - 1) != UNSURR)
317  dzp = 0 ;
318  p_dsdy->set_dzpuis(dzp) ;
319 
320  return *p_dsdy ;
321 
322 }
323 
324  //-----------------//
325  // d/dz //
326  //-----------------//
327 
328 const Scalar& Scalar::dsdz() const {
329 
330  // Protection
331  assert(etat != ETATNONDEF) ;
332 
333  // If the derivative has not been previously computed, the
334  // computation must be done by the appropriate routine of the mapping :
335 
336  if (p_dsdz == 0x0) {
337  p_dsdz = new Scalar(*mp) ;
338  if (etat == ETATUN) {
339  p_dsdz->set_etat_zero() ;
340  }
341  else {
343  }
344  }
345 
346  int dzp = (dzpuis == 0) ? 2 : dzpuis+1 ;
347  if (mp->get_mg()->get_type_r(mp->get_mg()->get_nzone() - 1) != UNSURR)
348  dzp = 0 ;
349  p_dsdz->set_dzpuis(dzp) ;
350 
351  return *p_dsdz ;
352 
353 }
354 
355  //-----------------//
356  // d/dx^i //
357  //-----------------//
358 
359 const Scalar& Scalar::deriv(int i) const {
360 
361  switch (i) {
362 
363  case 1 : {
364  return dsdx() ;
365  }
366 
367  case 2 : {
368  return dsdy() ;
369  }
370 
371  case 3 : {
372  return dsdz() ;
373  }
374 
375  default : {
376  cout << "Scalar::deriv : index i out of range !" << endl ;
377  cout << " i = " << i << endl ;
378  abort() ;
379  return dsdx() ; // Pour satisfaire le compilateur !
380  }
381 
382  }
383 
384 }
385 
386  //--------------------------//
387  // Covariant derivatives //
388  //--------------------------//
389 
390 const Vector& Scalar::derive_cov(const Metric& gam) const {
391 
392  const Vector* p_resu =
393  dynamic_cast<const Vector*>( &(Tensor::derive_cov(gam)) ) ;
394 
395  assert(p_resu != 0x0) ;
396 
397  return *p_resu ;
398 
399 }
400 
401 
402 const Vector& Scalar::derive_con(const Metric& gam) const {
403 
404  const Vector* p_resu =
405  dynamic_cast<const Vector*>( &(Tensor::derive_con(gam)) ) ;
406 
407  assert(p_resu != 0x0) ;
408 
409  return *p_resu ;
410 
411 }
412 
413 
414  //--------------------------//
415  // Lie derivative //
416  //--------------------------//
417 
418 
419 Scalar Scalar::derive_lie(const Vector& vv) const {
420 
421  Scalar resu(*mp) ;
422 
423  compute_derive_lie(vv, resu) ;
424 
425  return resu ;
426 
427 }
428 
429 
430 
431 
432  //---------------------//
433  // Laplacian //
434  //---------------------//
435 
436 const Scalar& Scalar::laplacian(int zec_mult_r) const {
437 
438  // Protection
439  assert(etat != ETATNONDEF) ;
440 
441  // If the Laplacian has not been previously computed, the
442  // computation must be done by the appropriate routine of the mapping :
443  if ( (p_lap == 0x0) || (zec_mult_r != ind_lap) ) {
444  if (p_lap != 0x0) {
445  delete p_lap ; // the Laplacian had been computed but with
446  // a different value of zec_mult_r
447  }
448  p_lap = new Scalar(*mp) ;
449  mp->laplacien(*this, zec_mult_r, *p_lap) ;
450  ind_lap = zec_mult_r ;
451  }
452 
453  return *p_lap ;
454 
455 }
456 
457  //-----------------------------//
458  // Angular Laplacian //
459  //-----------------------------//
460 
461 const Scalar& Scalar::lapang() const {
462 
463  // Protection
464  assert(etat != ETATNONDEF) ;
465 
466  // If the Laplacian has not been previously computed, the
467  // computation must be done by the appropriate routine of the mapping :
468  if ( p_lapang == 0x0 ) {
469  if (etat == ETATUN) {
470  p_lapang = new Scalar(*mp) ;
472  }
473  else {
474  p_lapang = new Scalar(*mp) ;
475  mp->lapang(*this, *p_lapang) ;
476  }
477  }
478 
480 
481  return *p_lapang ;
482 
483 }
484 
485 
486 
487  //---------------------//
488  // d/dradial //
489  //---------------------//
490 
491 const Scalar& Scalar::dsdradial() const {
492 
493  // Protection
494  assert(etat != ETATNONDEF) ;
495 
496  // If the derivative has not been previously computed, the
497  // computation must be done by the appropriate routine of the mapping :
498 
499  if (p_dsdradial == 0x0) {
500  p_dsdradial = new Scalar(*mp) ;
501  if (etat == ETATUN) {
503  }
504  else {
505  mp->dsdradial(*this, *p_dsdradial) ;
506  }
507  }
508 
509  int dzp = (dzpuis == 0) ? 2 : dzpuis+1 ;
510  if (mp->get_mg()->get_type_r(mp->get_mg()->get_nzone() - 1) != UNSURR)
511  dzp = 0 ;
512  p_dsdradial->set_dzpuis(dzp) ;
513 
514  return *p_dsdradial ;
515 
516 }
517 
518  //-----------------//
519  // d/drho //
520  //-----------------//
521 
522 const Scalar& Scalar::dsdrho() const {
523 
524  // Protection
525  assert(etat != ETATNONDEF) ;
526 
527  // If the derivative has not been previously computed, the
528  // computation must be done by the appropriate routine of the mapping :
529 
530  if (p_dsdrho == 0x0) {
531  p_dsdrho = new Scalar(*mp) ;
532  if (etat == ETATUN) {
534  }
535  else {
536  Scalar der_r (dsdr()) ;
537  Scalar der_t (srdsdt()) ;
538  Valeur val (der_r.get_spectral_va().mult_st() +
539  der_t.get_spectral_va().mult_ct()) ;
540 
541  Scalar res (*mp) ;
542  res = val ;
543 
544  *p_dsdrho = res ;
545  }
546  }
547 
548  int dzp = (dzpuis == 0) ? 2 : dzpuis+1 ;
549  if (mp->get_mg()->get_type_r(mp->get_mg()->get_nzone() - 1) != UNSURR)
550  dzp = 0 ;
551  p_dsdrho->set_dzpuis(dzp) ;
552 
553  return *p_dsdrho ;
554 
555 }
556 }
Scalar * p_srdsdt
Pointer on of *this (0x0 if not up to date)
Definition: scalar.h:422
Metric for tensor calculation.
Definition: metric.h:90
Scalar * p_srstdsdp
Pointer on of *this (0x0 if not up to date)
Definition: scalar.h:427
const Scalar & dsdradial() const
Returns of *this if the mapping is affine (class Map_af) and of *this if the mapping is logarithmic...
Definition: scalar_deriv.C:491
const Scalar & lapang() const
Returns the angular Laplacian of *this , where .
Definition: scalar_deriv.C:461
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
virtual void dsdt(const Scalar &uu, Scalar &resu) const =0
Computes of a Scalar .
Scalar(const Map &mpi)
Constructor from mapping.
Definition: scalar.C:210
Lorene prototypes.
Definition: app_hor.h:67
const Scalar & dsdt() const
Returns of *this .
Definition: scalar_deriv.C:208
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:777
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
Definition: scalar_deriv.C:436
virtual void dsdr(const Cmp &ci, Cmp &resu) const =0
Computes of a Cmp .
const Scalar & dsdz() const
Returns of *this , where .
Definition: scalar_deriv.C:328
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Scalar * p_dsdt
Pointer on of *this (0x0 if not up to date)
Definition: scalar.h:430
int ind_lap
Power of r by which the last computed Laplacian has been multiplied in the compactified external doma...
Definition: scalar.h:469
const Scalar & dsdx() const
Returns of *this , where .
Definition: scalar_deriv.C:266
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
int dzpuis
Power of r by which the quantity represented by this must be divided in the compactified external do...
Definition: scalar.h:409
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition: scalar_deriv.C:402
Tensor field of valence 1.
Definition: vector.h:188
virtual void comp_x_from_spherical(const Scalar &v_r, const Scalar &v_theta, const Scalar &v_phi, Scalar &v_x) const =0
Computes the Cartesian x component (with respect to bvect_cart ) of a vector given by its spherical c...
virtual void srdsdt(const Cmp &ci, Cmp &resu) const =0
Computes of a Cmp .
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:814
Scalar * p_lap
Pointer on the Laplacian of *this (0x0 if not up to date)
Definition: scalar.h:454
const Scalar & dsdrho() const
Returns of *this .
Definition: scalar_deriv.C:522
Scalar * p_dsdradial
Pointer on of *this.
Definition: scalar.h:461
virtual void comp_z_from_spherical(const Scalar &v_r, const Scalar &v_theta, Scalar &v_z) const =0
Computes the Cartesian z component (with respect to bvect_cart ) of a vector given by its spherical c...
void compute_derive_lie(const Vector &v, Tensor &resu) const
Computes the Lie derivative of this with respect to some vector field v (protected method; the public...
Scalar * p_dsdrho
Pointer on of *this.
Definition: scalar.h:464
const Scalar & dsdy() const
Returns of *this , where .
Definition: scalar_deriv.C:297
Scalar * p_dsdy
Pointer on of *this , where (0x0 if not up to date)
Definition: scalar.h:445
Scalar * p_lapang
Pointer on the Laplacian of *this (0x0 if not up to date)
Definition: scalar.h:458
virtual void srstdsdp(const Cmp &ci, Cmp &resu) const =0
Computes of a Cmp .
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
const Valeur & mult_st() const
Returns applied to *this.
const Scalar & deriv(int i) const
Returns of *this , where .
Definition: scalar_deriv.C:359
const Scalar & stdsdp() const
Returns of *this .
Definition: scalar_deriv.C:238
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor.C:1011
Scalar derive_lie(const Vector &v) const
Computes the derivative of this along a vector field v.
Definition: scalar_deriv.C:419
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition: tensor.C:1023
Scalar * p_dsdz
Pointer on of *this , where (0x0 if not up to date)
Definition: scalar.h:450
const Scalar & srdsdt() const
Returns of *this .
Definition: scalar_deriv.C:145
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary)...
Definition: scalar.h:402
Scalar * p_stdsdp
Pointer on of *this (0x0 if not up to date)
Definition: scalar.h:435
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
const Valeur & mult_ct() const
Returns applied to *this.
Scalar * p_dsdr
Pointer on of *this (0x0 if not up to date)
Definition: scalar.h:417
virtual void dsdradial(const Scalar &uu, Scalar &resu) const =0
Computes of a Scalar if the description is affine and if it is logarithmic.
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:491
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:301
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
virtual void comp_y_from_spherical(const Scalar &v_r, const Scalar &v_theta, const Scalar &v_phi, Scalar &v_y) const =0
Computes the Cartesian y component (with respect to bvect_cart ) of a vector given by its spherical c...
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607
const Scalar & srstdsdp() const
Returns of *this .
Definition: scalar_deriv.C:177
virtual void laplacien(const Scalar &uu, int zec_mult_r, Scalar &lap) const =0
Computes the Laplacian of a scalar field.
virtual void lapang(const Scalar &uu, Scalar &lap) const =0
Computes the angular Laplacian of a scalar field.
Scalar * p_dsdx
Pointer on of *this , where (0x0 if not up to date)
Definition: scalar.h:440
virtual void stdsdp(const Scalar &uu, Scalar &resu) const =0
Computes of a Scalar .