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