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
30char scalar_deriv_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_deriv.C,v 1.20 2014/10/13 08:53:46 j_novak Exp $" ;
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.20 2014/10/13 08:53:46 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
112namespace Lorene {
113const 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) {
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 ;
136
137 return *p_dsdr ;
138
139}
140
141 //--------------------//
142 // 1/r d/dtheta //
143 //--------------------//
144
145const 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 ;
167
168 return *p_srdsdt ;
169
170}
171
172
173 //------------------------------//
174 // 1/(r sin(theta)) d/dphi //
175 //------------------------------//
176
177const 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 ;
199
200 return *p_srstdsdp ;
201
202}
203
204 //--------------------//
205 // d/dtheta //
206 //--------------------//
207
208const 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) {
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
238const 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) {
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
266const 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) {
278 }
279 else {
280 mp->comp_x_from_spherical(dsdr(), srdsdt(), srstdsdp(), *p_dsdx) ;
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 ;
288
289 return *p_dsdx ;
290
291}
292
293 //-----------------//
294 // d/dy //
295 //-----------------//
296
297const 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) {
309 }
310 else {
311 mp->comp_y_from_spherical(dsdr(), srdsdt(), srstdsdp(), *p_dsdy) ;
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 ;
319
320 return *p_dsdy ;
321
322}
323
324 //-----------------//
325 // d/dz //
326 //-----------------//
327
328const 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) {
340 }
341 else {
342 mp->comp_z_from_spherical(dsdr(), srdsdt(), *p_dsdz) ;
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 ;
350
351 return *p_dsdz ;
352
353}
354
355 //-----------------//
356 // d/dx^i //
357 //-----------------//
358
359const 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
390const 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
402const 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
420
421 Scalar resu(*mp) ;
422
424
425 return resu ;
426
427}
428
429
430
431
432 //---------------------//
433 // Laplacian //
434 //---------------------//
435
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) ;
451 }
452
453 return *p_lap ;
454
455}
456
457 //-----------------------------//
458 // Angular Laplacian //
459 //-----------------------------//
460
461const 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
491const 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 ;
513
514 return *p_dsdradial ;
515
516}
517
518 //-----------------//
519 // d/drho //
520 //-----------------//
521
522const 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 ;
552
553 return *p_dsdrho ;
554
555}
556}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Metric for tensor calculation.
Definition metric.h:90
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Scalar & lapang() const
Returns the angular Laplacian of *this , where .
const Scalar & srdsdt() const
Returns of *this .
Scalar * p_lap
Pointer on the Laplacian of *this (0x0 if not up to date)
Definition scalar.h:448
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Scalar * p_dsdy
Pointer on of *this , where (0x0 if not up to date)
Definition scalar.h:439
const Scalar & deriv(int i) const
Returns of *this , where .
const Scalar & dsdz() const
Returns of *this , where .
const Scalar & dsdt() const
Returns of *this .
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
const Scalar & dsdy() const
Returns of *this , where .
Scalar * p_dsdrho
Pointer on of *this
Definition scalar.h:458
const Scalar & srstdsdp() const
Returns of *this .
Scalar * p_dsdr
Pointer on of *this (0x0 if not up to date)
Definition scalar.h:411
Scalar * p_stdsdp
Pointer on of *this (0x0 if not up to date)
Definition scalar.h:429
Scalar * p_dsdx
Pointer on of *this , where (0x0 if not up to date)
Definition scalar.h:434
const Scalar & dsdrho() const
Returns of *this .
const Scalar & stdsdp() const
Returns of *this .
Scalar * p_dsdz
Pointer on of *this , where (0x0 if not up to date)
Definition scalar.h:444
const Scalar & dsdx() const
Returns of *this , where .
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
const Scalar & dsdr() const
Returns of *this .
const Scalar & dsdradial() const
Returns of *this if the mapping is affine (class Map_af) and of *this if the mapping is logarithmic...
int ind_lap
Power of r by which the last computed Laplacian has been multiplied in the compactified external doma...
Definition scalar.h:463
Scalar * p_srstdsdp
Pointer on of *this (0x0 if not up to date)
Definition scalar.h:421
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary).
Definition scalar.h:396
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
Scalar * p_lapang
Pointer on the Laplacian of *this (0x0 if not up to date)
Definition scalar.h:452
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Scalar * p_dsdradial
Pointer on of *this
Definition scalar.h:455
Scalar * p_dsdt
Pointer on of *this (0x0 if not up to date)
Definition scalar.h:424
Scalar * p_srdsdt
Pointer on of *this (0x0 if not up to date)
Definition scalar.h:416
Scalar derive_lie(const Vector &v) const
Computes the derivative of this along a vector field v.
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
Tensor field of valence 1.
Definition vector.h:188
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
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:1014
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...
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition tensor.C:1002
Lorene prototypes.
Definition app_hor.h:64