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 char 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 
203 namespace Lorene {
204 Scalar::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 // -------------------------
215 Scalar::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 // ----------------
228 Scalar::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 //----------------------
238 Scalar::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 // ---------
248 Scalar::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() ;
282  va.set_etat_nondef() ;
284 
285 }
286 
287 void 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 
306 void Scalar::set_der_0x0() const {
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 
391 void 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 
436 void Scalar::operator=(const Tensor& uu) {
437 
438  assert(uu.valence == 0) ;
439 
440  operator=(*(uu.cmp[0])) ;
441 
442 }
443 
444 // From Scalar
445 // ----------
446 void Scalar::operator=(const Scalar& ci) {
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 // --------
499 void Scalar::operator=(const Cmp& ci) {
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 // -----------
549 void Scalar::operator=(const Valeur& vi) {
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 // ---------
592 void Scalar::operator=(const Mtbl& mi) {
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 // -----------
632 void 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 // --------
653 void Scalar::operator=(int n) {
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 //----------------------
674 Scalar::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 
686 void 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 // -----------
702 ostream& operator<<(ostream& ost, const Scalar& ci) {
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 
741 void Scalar::spectral_display(const char* comment,
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 
792 
793  va.std_base_scal_odd() ;
794 
795 }
796 
798 
799  va.set_base(bi) ;
800 
801 }
802 
803 
804  //--------------------------//
805  // dzpuis manipulations //
806  //--------------------------//
807 
808 void Scalar::set_dzpuis(int dzi) {
809 
810  dzpuis = dzi ;
811 
812 }
813 
814 bool 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 
873 bool 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 
890 double 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
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
Valeur va
The numerical value of the Cmp
Definition: cmp.h:464
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
int get_dzpuis() const
Returns dzpuis.
Definition: cmp.h:903
void set_dzpuis(int)
Set a value to dzpuis.
Definition: cmp.C:654
Base class for coordinate mappings.
Definition: map.h:670
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
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...
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
Base_val base
Bases of the spectral expansions.
Definition: mtbl_cf.h:200
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
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Definition: mtbl.h:132
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:299
int get_nzone() const
Gives the number of zones (domains)
Definition: mtbl.h:280
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
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:557
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
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:372
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: tbl.C:347
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
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
int get_etat() const
Returns the logical state.
Definition: valeur.h:726
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
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
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 coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Definition: valeur.h:480
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
int valence
Valence of the tensor (0 = scalar, 1 = vector, etc...)
Definition: tensor.h:298
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