LORENE
scalar.C
1/*
2 * Methods of class Scalar
3 *
4 * (see file scalar.h for documentation)
5 *
6 */
7
8/*
9 * Copyright (c) 2003-2004 Eric Gourgoulhon & Jerome Novak
10 *
11 * Copyright (c) 1999-2001 Philippe Grandclement (for preceding class Cmp)
12 *
13 * This file is part of LORENE.
14 *
15 * LORENE is free software; you can redistribute it and/or modify
16 * it under the terms of the GNU General Public License as published by
17 * the Free Software Foundation; either version 2 of the License, or
18 * (at your option) any later version.
19 *
20 * LORENE is distributed in the hope that it will be useful,
21 * but WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 * GNU General Public License for more details.
24 *
25 * You should have received a copy of the GNU General Public License
26 * along with LORENE; if not, write to the Free Software
27 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28 *
29 */
30
31
32char scalar_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar.C,v 1.40 2015/12/18 15:52:52 j_novak Exp $" ;
33
34
35/*
36 * $Id: scalar.C,v 1.40 2015/12/18 15:52:52 j_novak Exp $
37 * $Log: scalar.C,v $
38 * Revision 1.40 2015/12/18 15:52:52 j_novak
39 * New method is_nan() for class Scalar.
40 *
41 * Revision 1.39 2014/10/13 08:53:45 j_novak
42 * Lorene classes and functions now belong to the namespace Lorene.
43 *
44 * Revision 1.38 2014/10/06 15:16:15 j_novak
45 * Modified #include directives to use c++ syntax.
46 *
47 * Revision 1.37 2013/06/05 15:06:11 j_novak
48 * Legendre bases are treated as standard bases, when the multi-grid
49 * (Mg3d) is built with BASE_LEG.
50 *
51 * Revision 1.36 2013/03/27 09:51:42 e_gourgoulhon
52 * Restored log
53 *
54 *
55 * revision 1.35
56 * date: 2013/01/11 15:44:54; author: j_novak; state: Exp; lines: +15 -4
57 * Addition of Legendre bases (part 2).
58 *
59 * revision 1.34
60 * date: 2012/01/17 15:10:12; author: j_penner; state: Exp; lines: +2 -7
61 *
62 * revision 1.33
63 * date: 2012/01/17 10:26:48; author: j_penner; state: Exp; lines: +11 -3
64 * added a derivative with respect to the computational coordinate xi
65 *
66 * Revision 1.32 2011/04/08 12:55:50 e_gourgoulhon
67 * Scalar::val_point : division by r^dzpuis to return the actual
68 * (i.e. physical) value of the field in the external compactified
69 * domain.
70 *
71 * Revision 1.31 2005/10/25 08:56:39 p_grandclement
72 * addition of std_spectral_base in the case of odd functions near the origin
73 *
74 * Revision 1.30 2005/03/11 13:16:37 j_novak
75 * Corrected an error in multipole_spectrum().
76 *
77 * Revision 1.29 2004/10/11 15:09:04 j_novak
78 * The radial manipulation functions take Scalar as arguments, instead of Cmp.
79 * Added a conversion operator from Scalar to Cmp.
80 * The Cmp radial manipulation function make conversion to Scalar, call to the
81 * Map_radial version with a Scalar argument and back.
82 *
83 * Revision 1.28 2004/08/24 09:14:51 p_grandclement
84 * Addition of some new operators, like Poisson in 2d... It now requieres the
85 * GSL library to work.
86 *
87 * Also, the way a variable change is stored by a Param_elliptic is changed and
88 * no longer uses Change_var but rather 2 Scalars. The codes using that feature
89 * will requiere some modification. (It should concern only the ones about monopoles)
90 *
91 * Revision 1.27 2004/06/22 08:50:00 p_grandclement
92 * Addition of everything needed for using the logarithmic mapping
93 *
94 * Revision 1.26 2004/06/11 14:29:58 j_novak
95 * Scalar::multipole_spectrum() is now a const method.
96 *
97 * Revision 1.25 2004/03/04 09:51:36 e_gourgoulhon
98 * Method dz_nonzero() : treated the case ETATUN.
99 *
100 * Revision 1.24 2004/02/19 22:11:50 e_gourgoulhon
101 * Added argument "comment" in method spectral_display.
102 *
103 * Revision 1.23 2004/01/12 15:32:25 j_novak
104 * Yet another problem with ETATUN
105 *
106 * Revision 1.22 2003/11/05 15:31:13 e_gourgoulhon
107 * Method set_etat_qcq(): del_t() is not invoqued for etat == ETATUN.
108 *
109 * Revision 1.21 2003/11/03 10:25:27 j_novak
110 * Suppressed the call to del_deriv() in set_etat_qcq() method.
111 *
112 * Revision 1.20 2003/10/28 21:33:13 e_gourgoulhon
113 * In operator=(const Scalar& ci) changed the place of va.del_t()
114 * in order to correct an error in the case where both this and ci
115 * are in the state ETATUN.
116 *
117 * Revision 1.19 2003/10/28 12:36:10 e_gourgoulhon
118 * Corrected operator=(const Scalar& ci) for the case ETATUN: the
119 * spectral bases are now copied.
120 *
121 * Revision 1.18 2003/10/27 14:51:25 e_gourgoulhon
122 * Correction of errors in constructor from a Tensor.
123 *
124 * Revision 1.17 2003/10/22 08:35:30 j_novak
125 * Error corrected
126 *
127 * Revision 1.16 2003/10/19 19:54:37 e_gourgoulhon
128 * -- Modified method spectral_display: now calling Valeur::display_coef.
129 *
130 * Revision 1.15 2003/10/15 16:03:38 j_novak
131 * Added the angular Laplace operator for Scalar.
132 *
133 * Revision 1.14 2003/10/15 10:43:06 e_gourgoulhon
134 * Added new members p_dsdt and p_stdsdp.
135 *
136 * Revision 1.13 2003/10/13 13:52:40 j_novak
137 * Better managment of derived quantities.
138 *
139 * Revision 1.12 2003/10/11 14:42:16 e_gourgoulhon
140 * Suppressed unusued argument new_triad in method change_triad.
141 *
142 * Revision 1.11 2003/10/10 15:57:29 j_novak
143 * Added the state one (ETATUN) to the class Scalar
144 *
145 * Revision 1.10 2003/10/07 08:05:03 j_novak
146 * Added an assert for the constructor from a Tensor.
147 *
148 * Revision 1.9 2003/10/06 16:16:03 j_novak
149 * New constructor from a Tensor.
150 *
151 * Revision 1.8 2003/10/06 13:58:48 j_novak
152 * The memory management has been improved.
153 * Implementation of the covariant derivative with respect to the exact Tensor
154 * type.
155 *
156 * Revision 1.7 2003/10/05 21:15:42 e_gourgoulhon
157 * - Suppressed method std_spectral_base_scal().
158 * - Added method std_spectral_base().
159 *
160 * Revision 1.6 2003/10/01 13:04:44 e_gourgoulhon
161 * The method Tensor::get_mp() returns now a reference (and not
162 * a pointer) onto a mapping.
163 *
164 * Revision 1.5 2003/09/29 12:52:58 j_novak
165 * Methods for changing the triad are implemented.
166 *
167 * Revision 1.4 2003/09/24 20:55:27 e_gourgoulhon
168 * Added -- constructor by conversion of a Cmp
169 * -- operator=(const Cmp&)
170 *
171 * Revision 1.3 2003/09/24 15:10:54 j_novak
172 * Suppression of the etat flag in class Tensor (still present in Scalar)
173 *
174 * Revision 1.2 2003/09/24 10:21:07 e_gourgoulhon
175 * added more methods
176 *
177 * Revision 1.1 2003/09/23 08:52:17 e_gourgoulhon
178 * new version
179 *
180 *
181 * $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar.C,v 1.40 2015/12/18 15:52:52 j_novak Exp $
182 *
183 */
184
185// headers C
186#include <cassert>
187#include <cstdlib>
188#include <cmath>
189
190// headers Lorene
191#include "tensor.h"
192#include "type_parite.h"
193#include "utilitaires.h"
194#include "proto.h"
195#include "cmp.h"
196
197
198 //---------------//
199 // Constructors //
200 //---------------//
201
202
203namespace Lorene {
204Scalar::Scalar(const Map& mpi) : Tensor(mpi), etat(ETATNONDEF), dzpuis(0),
205 va(mpi.get_mg()) {
206
207 cmp[0] = this ;
208 set_der_0x0() ;
209
210}
211
212
213// Constructor from a Tensor
214// -------------------------
215Scalar::Scalar(const Tensor& ti) : Tensor(*(ti.mp)), etat(ti.cmp[0]->etat),
216 dzpuis(ti.cmp[0]->dzpuis), va(ti.cmp[0]->va) {
217
218 assert(ti.valence == 0) ;
219
220 cmp[0] = this ;
221 set_der_0x0() ;
222
223}
224
225
226// Copy constructor
227// ----------------
228Scalar::Scalar(const Scalar& sci) : Tensor(*(sci.mp)), etat(sci.etat),
229 dzpuis(sci.dzpuis), va(sci.va) {
230
231 cmp[0] = this ;
232 set_der_0x0() ; // On ne recopie pas les derivees
233
234}
235
236// Conversion from a Cmp
237//----------------------
238Scalar::Scalar(const Cmp& ci) : Tensor(*(ci.get_mp())),
239 etat(ci.get_etat()),
240 dzpuis(ci.get_dzpuis()),
241 va(ci.va) {
242 cmp[0] = this ;
243 set_der_0x0() ;
244}
245
246// From file
247// ---------
248Scalar::Scalar(const Map& mpi, const Mg3d& mgi, FILE* fd) : Tensor(mpi),
249 va(mgi, fd) {
250
251 assert( mpi.get_mg() == &mgi ) ;
252
253 fread_be(&etat, sizeof(int), 1, fd) ; // L'etat
254 fread_be(&dzpuis, sizeof(int), 1, fd) ; // dzpuis
255
256 cmp[0] = this ;
257
258 set_der_0x0() ; // Les derivees sont initialisees a zero
259
260}
261
262 //--------------//
263 // Destructor //
264 //--------------//
265
266// Destructeur
268 del_t() ;
269 cmp[0] = 0x0 ; //cmp[0] was set to 'this' before (now 0x0 to break a
270 // in loop in ~Tensor!)
271
272}
273
274 //-----------------------//
275 // Memory management //
276 //-----------------------//
277
278// Destructeur logique
280
281 va.del_t() ;
284
285}
286
287void Scalar::del_deriv() const{
288 if (p_dsdr != 0x0) delete p_dsdr ;
289 if (p_srdsdt != 0x0) delete p_srdsdt ;
290 if (p_srstdsdp != 0x0) delete p_srstdsdp ;
291 if (p_dsdt != 0x0) delete p_dsdt ;
292 if (p_stdsdp != 0x0) delete p_stdsdp ;
293 if (p_dsdx != 0x0) delete p_dsdx ;
294 if (p_dsdy != 0x0) delete p_dsdy ;
295 if (p_dsdz != 0x0) delete p_dsdz ;
296 if (p_lap != 0x0) delete p_lap ;
297 if (p_lapang != 0x0) delete p_lapang ;
298 if (p_integ != 0x0) delete p_integ ;
299 if (p_dsdradial != 0x0) delete p_dsdradial ;
300 if (p_dsdrho != 0x0) delete p_dsdrho ;
301 set_der_0x0() ;
302
304}
305
307 p_dsdr = 0x0 ;
308 p_srdsdt = 0x0 ;
309 p_srstdsdp = 0x0 ;
310 p_dsdt = 0x0 ;
311 p_stdsdp = 0x0 ;
312 p_dsdx = 0x0 ;
313 p_dsdy = 0x0 ;
314 p_dsdz = 0x0 ;
315 p_lap = 0x0 ;
316 p_lapang = 0x0 ;
317 ind_lap = - 1 ;
318 p_integ = 0x0 ;
319 p_dsdradial = 0x0 ;
320 p_dsdrho = 0x0 ;
321}
322
323// ETATZERO
325 if (etat == ETATZERO) return ;
326 else {
327 del_deriv() ;
328 va.set_etat_zero() ;
329 etat = ETATZERO ;
330 }
331}
332
333// ETATUN
335 if (etat == ETATUN) return ;
336 else {
337 del_deriv() ;
338 va = 1 ;
339 etat = ETATUN ;
340 }
341}
342
343// ETATNONDEF
345 if (etat == ETATNONDEF) return ;
346 else {
347 del_t() ;
348 etat = ETATNONDEF ;
349 }
350}
351
352// ETATQCQ
354
355 if (etat == ETATQCQ) return ;
356 else {
357 if (etat != ETATUN) del_t() ;
358 etat = ETATQCQ ;
359 }
360}
361
362
363// Allocates everything
364// --------------------
366
367 set_etat_qcq() ;
368 va.set_etat_c_qcq() ; // allocation in configuration space
369 Mtbl* mt = va.c ;
370 mt->set_etat_qcq() ;
371 for (int l=0; l<mt->get_nzone(); l++) {
372 mt->t[l]->set_etat_qcq() ;
373 }
374
375}
376
377
378
379// ZERO hard
381
382 va.annule_hard() ;
383 del_deriv() ;
384 etat = ETATQCQ ;
385}
386
387
388// Sets the Scalar to zero in several domains
389// ---------------------------------------
390
391void Scalar::annule(int l_min, int l_max) {
392
393 // Cas particulier: annulation globale :
394 if ( (l_min == 0) && (l_max == va.mg->get_nzone()-1) ) {
395 set_etat_zero() ;
396 return ;
397 }
398
399 assert( etat != ETATNONDEF ) ;
400
401 if ( etat == ETATZERO ) {
402 return ; // rien n'a faire si c'est deja zero
403 }
404 else {
405 assert( (etat == ETATQCQ) || (etat == ETATUN) ) ; // sinon...
406
407 etat = ETATQCQ ;
408 va.annule(l_min, l_max) ; // Annule la Valeur
409
410 // Annulation des membres derives
411 if (p_dsdr != 0x0) p_dsdr->annule(l_min, l_max) ;
412 if (p_srdsdt != 0x0) p_srdsdt->annule(l_min, l_max) ;
413 if (p_srstdsdp != 0x0) p_srstdsdp->annule(l_min, l_max) ;
414 if (p_dsdt != 0x0) p_dsdt->annule(l_min, l_max) ;
415 if (p_stdsdp != 0x0) p_stdsdp->annule(l_min, l_max) ;
416 if (p_dsdx != 0x0) p_dsdx->annule(l_min, l_max) ;
417 if (p_dsdy != 0x0) p_dsdy->annule(l_min, l_max) ;
418 if (p_dsdz != 0x0) p_dsdz->annule(l_min, l_max) ;
419 if (p_lap != 0x0) p_lap->annule(l_min, l_max) ;
420 if (p_lapang != 0x0) p_lapang->annule(l_min, l_max) ;
421 if (p_integ != 0x0) delete p_integ ;
422 if (p_dsdradial != 0x0) p_dsdradial->annule(l_min, l_max) ;
423 }
424
425}
426
427
428 //------------//
429 // Assignment //
430 //------------//
431
432
433// From tensor
434// -----------
435
437
438 assert(uu.valence == 0) ;
439
440 operator=(*(uu.cmp[0])) ;
441
442}
443
444// From Scalar
445// ----------
447
448
449 assert(&ci != this) ; // pour eviter l'auto-affectation
450
451 // Les elements fixes
452 assert( mp == ci.mp ) ;
453
454 dzpuis = ci.dzpuis ;
455
456 // La valeur eventuelle
457 switch(ci.etat) {
458 case ETATNONDEF: {
459 set_etat_nondef() ;
460 break ; // valeur par defaut
461 }
462
463 case ETATZERO: {
464 set_etat_zero() ;
465 break ;
466 }
467
468 case ETATUN: {
469 set_etat_one() ;
470 va.set_base( ci.va.get_base() ) ;
471 break ;
472 }
473
474 case ETATQCQ: {
475 // Menage general de la Valeur, mais pas des quantites derivees !
476 va.del_t() ;
477
478 set_etat_qcq() ; // set_etat_qcq n'appelle pas del_deriv()
479 va = ci.va ;
480
481 // On detruit les quantites derivees (seulement lorsque tout est fini !)
482 del_deriv() ;
483
484 break ;
485 }
486
487 default: {
488 cout << "Unkwown state in Scalar::operator=(const Scalar&) !"
489 << endl ;
490 abort() ;
491 break ;
492 }
493 }
494
495}
496
497// From Cmp
498// --------
500
501
502 // Menage general de la Valeur, mais pas des quantites derivees !
503 va.del_t() ;
504
505 // Les elements fixes
506 assert( mp == ci.get_mp() ) ;
507
508 dzpuis = ci.get_dzpuis() ;
509
510 // La valeur eventuelle
511 switch(ci.get_etat()) {
512
513 case ETATNONDEF: {
514 set_etat_nondef() ;
515 break ; // valeur par defaut
516 }
517
518 case ETATZERO: {
519 set_etat_zero() ;
520 break ;
521 }
522
523 case ETATQCQ: {
524 set_etat_qcq() ;
525 va = ci.va ;
526
527 // On detruit les quantites derivees (seulement lorsque tout est fini !)
528 del_deriv() ;
529
530 break ;
531 }
532
533 default: {
534 cout << "Unkwown state in Scalar::operator=(const Cmp&) !"
535 << endl ;
536 abort() ;
537 break ;
538 }
539 }
540
541}
542
543
544
545
546
547// From Valeur
548// -----------
550
551 // Traitement de l'auto-affectation :
552 if (&vi == &va) {
553 return ;
554 }
555
556 // Protection
557 assert(vi.get_etat() != ETATNONDEF) ;
558
559 // Menage general de la Valeur, mais pas des quantites derivees !
560 va.del_t() ;
561
562
563 // La valeure eventuelle
564 switch(vi.get_etat()) {
565
566 case ETATZERO: {
567 set_etat_zero() ;
568 break ;
569 }
570
571 case ETATQCQ: {
572 set_etat_qcq() ;
573 va = vi ;
574
575 // On detruit les quantites derivees (seulement lorsque tout est fini !)
576 del_deriv() ;
577
578 break ;
579 }
580
581 default: {
582 cout << "Unkwown state in Scalar::operator=(const Valeur&) !" << endl ;
583 abort() ;
584 break ;
585 }
586 }
587
588}
589
590// From Mtbl
591// ---------
593
594 // Protection
595 assert(mi.get_etat() != ETATNONDEF) ;
596
597 assert(&mi != va.c) ; // pour eviter l'auto-affectation
598
599
600 // Menage general de la Valeur, mais pas des quantites derivees !
601 va.del_t() ;
602
603 // La valeure eventuelle
604 switch(mi.get_etat()) {
605 case ETATZERO: {
606 set_etat_zero() ;
607 break ;
608 }
609
610 case ETATQCQ: {
611 set_etat_qcq() ;
612 va = mi ;
613
614 // On detruit les quantites derivees (seulement lorsque tout est fini !)
615 del_deriv() ;
616
617 break ;
618 }
619
620 default: {
621 cout << "Unkwown state in Scalar::operator=(const Mtbl&) !" << endl ;
622 abort() ;
623 break ;
624 }
625 }
626
627
628}
629
630// From double
631// -----------
632void Scalar::operator=(double x) {
633
634 if (x == double(0)) {
635 set_etat_zero() ;
636 }
637 else {
638 if (x == double(1)) {
639 set_etat_one() ;
640 }
641 else {
642 set_etat_qcq() ;
643 del_deriv() ;
644 va = x ;
645 }
646 }
647
648 dzpuis = 0 ;
649}
650
651// From int
652// --------
654
655 if (n == 0) {
656 set_etat_zero() ;
657 }
658 else {
659 if (n == 1) {
660 set_etat_one() ;
661 }
662 else {
663 set_etat_qcq() ;
664 del_deriv() ;
665 va = n ;
666 }
667 }
668 dzpuis = 0 ;
669
670}
671
672// Conversion to a Cmp
673//----------------------
674Scalar::operator Cmp() const {
675
676 Cmp resu(mp) ;
677 resu = va ;
678 resu.set_dzpuis(dzpuis) ;
679 return resu ;
680
681}
682 //------------//
683 // Sauvegarde //
684 //------------//
685
686void Scalar::sauve(FILE* fd) const {
687
688 va.sauve(fd) ; // la valeur (en premier pour la construction
689 // lors de la lecture du fichier)
690
691 fwrite_be(&etat, sizeof(int), 1, fd) ; // l'etat
692 fwrite_be(&dzpuis, sizeof(int), 1, fd) ; // dzpuis
693
694}
695
696 //------------//
697 // Impression //
698 //------------//
699
700// Operator <<
701// -----------
703
704 switch(ci.etat) {
705 case ETATNONDEF: {
706 ost << "*** UNDEFINED STATE *** " << endl ;
707 break ;
708 }
709
710 case ETATZERO: {
711 ost << "*** identically ZERO ***" << endl ;
712 break ;
713 }
714
715 case ETATUN: {
716 ost << "*** identically ONE ***" << endl ;
717 break ;
718 }
719
720 case ETATQCQ: {
721 ost << "*** dzpuis = " << ci.get_dzpuis() << endl ;
722 ost << ci.va << endl ;
723 break ;
724 }
725
726 default: {
727 cout << "operator<<(ostream&, const Scalar&) : unknown state !"
728 << endl ;
729 abort() ;
730 break ;
731 }
732 }
733
734 // Termine
735 return ost ;
736}
737
738// spectral_display
739//-----------------
740
742 double thres, int precis, ostream& ost) const {
743
744
745 if (comment != 0x0) {
746 ost << comment << " : " << endl ;
747 }
748
749 // Cas particuliers
750 //-----------------
751
752 if (etat == ETATNONDEF) {
753 ost << "*** UNDEFINED ***" << endl << endl ;
754 return ;
755 }
756
757 if (etat == ETATZERO) {
758 ost << "*** identically ZERO ***" << endl ;
759 return ;
760 }
761
762 if (etat == ETATUN) {
763 ost << "*** identically ONE ***" << endl ;
764 return ;
765 }
766
767 // Cas general : on affiche la Valeur
768 //------------
769
770 if (dzpuis != 0) {
771 ost << "*** dzpuis = " << dzpuis << endl ;
772 }
773 va.display_coef(thres, precis, ost) ;
774
775}
776
777
778
779
780 //------------------------------------//
781 // Spectral bases of the Valeur va //
782 //------------------------------------//
783
785
786 va.std_base_scal() ;
787
788}
789
790
796
798
799 va.set_base(bi) ;
800
801}
802
803
804 //--------------------------//
805 // dzpuis manipulations //
806 //--------------------------//
807
809
810 dzpuis = dzi ;
811
812}
813
814bool Scalar::dz_nonzero() const {
815
816 assert(etat != ETATNONDEF) ;
817
818 const Mg3d* mg = mp->get_mg() ;
819
820 int nzm1 = mg->get_nzone() - 1;
821 if (mg->get_type_r(nzm1) != UNSURR) {
822 return false ;
823 }
824
825 if (etat == ETATZERO) {
826 return false ;
827 }
828
829 if (etat == ETATUN) {
830 return true ;
831 }
832
833 assert( etat == ETATQCQ ) ;
834
835 if (va.etat == ETATZERO) {
836 return false ;
837 }
838
839 assert(va.etat == ETATQCQ) ;
840
841 if (va.c != 0x0) {
842 if ( (va.c)->get_etat() == ETATZERO ) {
843 return false ;
844 }
845
846 assert( (va.c)->get_etat() == ETATQCQ ) ;
847 if ( (va.c)->t[nzm1]->get_etat() == ETATZERO ) {
848 return false ;
849 }
850 else {
851 assert( (va.c)->t[nzm1]->get_etat() == ETATQCQ ) ;
852 return true ;
853 }
854 }
855 else{
856 assert(va.c_cf != 0x0) ;
857 if ( (va.c_cf)->get_etat() == ETATZERO ) {
858 return false ;
859 }
860 assert( (va.c_cf)->get_etat() == ETATQCQ ) ;
861 if ( (va.c_cf)->t[nzm1]->get_etat() == ETATZERO ) {
862 return false ;
863 }
864 else {
865 assert( (va.c_cf)->t[nzm1]->get_etat() == ETATQCQ ) ;
866 return true ;
867 }
868
869 }
870
871}
872
873bool Scalar::check_dzpuis(int dzi) const {
874
875 if (dz_nonzero()) { // the check must be done
876 return (dzpuis == dzi) ;
877 }
878 else{
879 return true ;
880 }
881
882}
883
884
885
886 //---------------------------------------//
887 // Value at an arbitrary point //
888 //---------------------------------------//
889
890double Scalar::val_point(double r, double theta, double phi) const {
891
892
893 assert(etat != ETATNONDEF) ;
894
895 if (etat == ETATZERO) {
896 return double(0) ;
897 }
898
899 if (etat == ETATUN) {
900 return double(1) ;
901 }
902
903 assert(etat == ETATQCQ) ;
904
905 // 1/ Search for the domain and the grid coordinates (xi,theta',phi')
906 // which corresponds to the point (r,theta,phi)
907
908 int l ;
909 double xi ;
910 mp->val_lx(r, theta, phi, l, xi) ; // call of val_lx with default
911 // accuracy parameters
912 // 2/ Call to the Valeur version
913 if (l == mp->get_mg()->get_nzone() - 1){ // in the last domain, one must take into account dzpuis
914 switch (dzpuis) {
915 case 0:
916 {
917 return va.val_point(l, xi, theta, phi);
918 break;
919 }
920 case 1:
921 {
922 return va.val_point(l, xi, theta, phi)/r ;
923 break;
924 }
925 case 2:
926 {
927 return va.val_point(l, xi, theta, phi)/(r*r) ;
928 break;
929 }
930 case 3:
931 {
932 return va.val_point(l, xi, theta, phi)/(r*r*r) ;
933 break;
934 }
935 case 4:
936 {
937 return va.val_point(l, xi, theta, phi)/(r*r*r*r) ;
938 break;
939 }
940 default:
941 {
942 cout << "scalar::val_point unexpected value of dzpuis !" << endl;
943 abort();
944 }
945 }
946 }
947 else{
948 return va.val_point(l, xi, theta, phi);
949 }
950}
951
952 //---------------------------------//
953 // Multipolar spectrum //
954 //---------------------------------//
955
957
958 assert (etat != ETATNONDEF) ;
959
960 const Mg3d* mg = mp->get_mg() ;
961 int nzone = mg->get_nzone() ;
962 int lmax = 0 ;
963
964 for (int lz=0; lz<nzone; lz++)
965 lmax = (lmax < 2*mg->get_nt(lz) - 1 ? 2*mg->get_nt(lz) - 1 : lmax) ;
966
967 Tbl resu(nzone, lmax) ;
968 if (etat == ETATZERO) {
969 resu.set_etat_zero() ;
970 return resu ;
971 }
972
973 assert((etat == ETATQCQ) || (etat == ETATUN));
974
975 Valeur va_copy = va ;
976
977 va_copy.coef() ;
978 va_copy.ylm() ;
979 resu.annule_hard() ;
980 const Base_val& base = va_copy.c_cf->base ;
981 int m_quant, l_quant, base_r ;
982 for (int lz=0; lz<nzone; lz++)
983 for (int k=0 ; k<mg->get_np(lz) ; k++)
984 for (int j=0 ; j<mg->get_nt(lz) ; j++) {
985 if (nullite_plm(j, mg->get_nt(lz), k, mg->get_np(lz), base) == 1)
986 {
987 // quantic numbers and spectral bases
988 donne_lm(nzone, lz, j, k, base, m_quant, l_quant, base_r) ;
989 for (int i=0; i<mg->get_nr(lz); i++) resu.set(lz, l_quant)
990 += fabs((*va_copy.c_cf)(lz, k, j, i)) ;
991 }
992 }
993
994 return resu ;
995}
996
998
999 cout << "WARNING: Scalar::change_triad : "<< endl ;
1000 cout << "This method does nothing ... " << endl ;
1001
1002}
1003
1004 //---------------------------------//
1005 // Check for NaNs //
1006 //---------------------------------//
1007
1008 bool Scalar::is_nan(bool verb) const {
1009
1010 bool resu = false ;
1011 if (etat == ETATQCQ) {
1012 bool phys = false ;
1013 bool coef = false ;
1014 if (va.c != 0x0)
1015 phys = true ;
1016 if (va.c_cf != 0x0)
1017 coef = true ;
1018 int nz = mp->get_mg()->get_nzone() ;
1019 for (int lz=0; lz<nz; lz++) {
1020 bool domain_p = false ;
1021 bool domain_c = false ;
1022 if (phys) domain_p = (va.c->get_etat() == ETATQCQ) ;
1023 if (coef) domain_c = (va.c_cf->get_etat() == ETATQCQ) ;
1024 if (domain_p) {
1025 int np = mp->get_mg()->get_np(lz) ;
1026 int nt = mp->get_mg()->get_nt(lz) ;
1027 int nr = mp->get_mg()->get_nr(lz) ;
1028 for (int k=0; k<np; k++) {
1029 for (int j=0; j<nt; j++) {
1030 for (int i=0; i<nr; i++) {
1031 if (isnan( (*va.c)(lz, k, j, i) ) ) {
1032 resu = true ;
1033 if (verb) {
1034 cout << "NaN found at physical grid point domain = " << lz
1035 << ", i= " << i << ", j= " << j << ", k= " << k << endl ;
1036 }
1037 }
1038 } // i loop
1039 } // j loop
1040 } // k loop
1041 }
1042 if (domain_c) {
1043 int np = mp->get_mg()->get_np(lz) + 2 ;
1044 int nt = mp->get_mg()->get_nt(lz) ;
1045 int nr = mp->get_mg()->get_nr(lz) ;
1046 for (int k=0; k<np; k++) {
1047 for (int j=0; j<nt; j++) {
1048 for (int i=0; i<nr; i++) {
1049 if (isnan( (*va.c_cf)(lz, k, j, i) ) ) {
1050 resu = true ;
1051 if (verb) {
1052 cout << "NaN found at coefficient, domain = " << lz
1053 << ", i= " << i << ", j= " << j << ", k= " << k << endl ;
1054 }
1055 }
1056 } // i loop
1057 } // j loop
1058 } // k loop
1059 }
1060 } // lz loop
1061 } // etat condition
1062
1063 return resu ;
1064
1065 }
1066
1067
1068}
Bases of the spectral expansions.
Definition base_val.h:322
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition base_vect.h:105
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Base class for coordinate mappings.
Definition map.h:670
Multi-domain grid.
Definition grilles.h:273
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:474
int get_etat() const
Returns the logical state.
Definition mtbl_cf.h:456
Multi-domain array.
Definition mtbl.h:118
int get_etat() const
Gives the logical state.
Definition mtbl.h:277
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
virtual void sauve(FILE *) const
Save in a file.
Definition scalar.C:686
void set_etat_one()
Sets the logical state to ETATUN (one).
Definition scalar.C:334
virtual void del_deriv() const
Logical destructor of the derivatives.
Definition scalar.C:287
Scalar * p_lap
Pointer on the Laplacian of *this (0x0 if not up to date)
Definition scalar.h:448
Scalar * p_dsdy
Pointer on of *this , where (0x0 if not up to date)
Definition scalar.h:439
virtual ~Scalar()
Destructor.
Definition scalar.C:267
Scalar * p_dsdrho
Pointer on of *this
Definition scalar.h:458
void operator=(const Scalar &a)
Assignment to another Scalar defined on the same mapping.
Definition scalar.C:446
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
bool dz_nonzero() const
Returns true if the last domain is compactified and *this is not zero in this domain.
Definition scalar.C:814
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
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:784
bool is_nan(bool verbose=false) const
Looks for NaNs (not a number) in the scalar field.
Definition scalar.C:1008
Scalar * p_dsdz
Pointer on of *this , where (0x0 if not up to date)
Definition scalar.h:444
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
virtual void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition scalar.C:997
int ind_lap
Power of r by which the last computed Laplacian has been multiplied in the compactified external doma...
Definition scalar.h:463
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition scalar.C:391
Scalar * p_srstdsdp
Pointer on of *this (0x0 if not up to date)
Definition scalar.h:421
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
void del_t()
Logical destructor.
Definition scalar.C:279
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary).
Definition scalar.h:396
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition scalar.C:365
virtual void set_etat_nondef()
Sets the logical state to ETATNONDEF (undefined).
Definition scalar.C:344
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
double val_point(double r, double theta, double phi) const
Computes the value of the field at an arbitrary point , by means of the spectral expansion.
Definition scalar.C:890
Tbl multipole_spectrum() const
Gives the spectrum in terms of multipolar modes l .
Definition scalar.C:956
virtual void std_spectral_base_odd()
Sets the spectral bases of the Valeur va to the standard odd ones for a scalar field.
Definition scalar.C:791
Valeur va
The numerical value of the Scalar
Definition scalar.h:405
virtual void spectral_display(const char *comment=0x0, double threshold=1.e-7, int precision=4, ostream &ostr=cout) const
Displays the spectral coefficients and the associated basis functions.
Definition scalar.C:741
Scalar * p_dsdradial
Pointer on of *this
Definition scalar.h:455
Tbl * p_integ
Pointer on the space integral of *this (values in each domain) (0x0 if not up to date)
Definition scalar.h:468
Scalar * p_dsdt
Pointer on of *this (0x0 if not up to date)
Definition scalar.h:424
void set_der_0x0() const
Sets the pointers for derivatives to 0x0.
Definition scalar.C:306
Scalar * p_srdsdt
Pointer on of *this (0x0 if not up to date)
Definition scalar.h:416
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
Definition scalar.C:797
int dzpuis
Power of r by which the quantity represented by this must be divided in the compactified external d...
Definition scalar.h:403
Basic array class.
Definition tbl.h:161
Tensor handling.
Definition tensor.h:288
Values and coefficients of a (real-value) function.
Definition valeur.h:287
void del_t()
Logical destructor.
Definition valeur.C:626
double val_point(int l, double x, double theta, double phi) const
Computes the value of the field represented by *this at an arbitrary point, by means of the spectral ...
Definition valeur.C:882
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c ).
Definition valeur.C:701
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition valeur.C:689
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition valeur.C:810
Mtbl * c
Values of the function at the points of the multi-grid
Definition valeur.h:299
void display_coef(double threshold=1.e-7, int precision=4, ostream &ostr=cout) const
Displays the spectral coefficients and the associated basis functions.
Definition valeur.C:536
void annule(int l)
Sets the Valeur to zero in a given domain.
Definition valeur.C:744
const Mg3d * mg
Multi-grid Mgd3 on which this is defined.
Definition valeur.h:292
void std_base_scal_odd()
Sets the bases for spectral expansions (member base ) to the standard odd ones for a scalar.
Definition valeur.C:830
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Definition valeur.C:824
void annule_hard()
Sets the Valeur to zero in a hard way.
Definition valeur.C:723
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition valeur.h:295
void set_etat_nondef()
Sets the logical state to ETATNONDEF (undefined).
Definition valeur.C:695
void sauve(FILE *) const
Save in a file.
Definition valeur.C:475
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition fread_be.C:69
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition fwrite_be.C:70
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition tensor.h:315
virtual void del_deriv() const
Deletes the derived quantities.
Definition tensor.C:398
Lorene prototypes.
Definition app_hor.h:64