LORENE
valeur_math.C
1 /*
2  * Mathematical functions for class Valeur
3  *
4  */
5 
6 /*
7  * Copyright (c) 1999-2000 Jean-Alain Marck
8  * Copyright (c) 1999-2001 Eric Gourgoulhon
9  * Copyright (c) 1999-2001 Philippe Grandclement
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 
32 /*
33  * $Id: valeur_math.C,v 1.6 2016/12/05 16:18:20 j_novak Exp $
34  * $Log: valeur_math.C,v $
35  * Revision 1.6 2016/12/05 16:18:20 j_novak
36  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
37  *
38  * Revision 1.5 2014/10/13 08:53:50 j_novak
39  * Lorene classes and functions now belong to the namespace Lorene.
40  *
41  * Revision 1.4 2014/10/06 15:13:23 j_novak
42  * Modified #include directives to use c++ syntax.
43  *
44  * Revision 1.3 2012/01/17 10:39:27 j_penner
45  * added a Heaviside function
46  *
47  * Revision 1.2 2002/10/16 14:37:15 j_novak
48  * Reorganization of #include instructions of standard C++, in order to
49  * use experimental version 3 of gcc.
50  *
51  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
52  * LORENE
53  *
54  * Revision 2.4 1999/12/02 17:57:42 phil
55  * *** empty log message ***
56  *
57  * Revision 2.3 1999/11/09 16:18:02 phil
58  * ajout de racine_cubique
59  *
60  * Revision 2.2 1999/10/29 15:14:50 eric
61  * Ajout de fonctions mathematiques (abs, norme, etc...).
62  *
63  * Revision 2.1 1999/03/01 15:11:03 eric
64  * *** empty log message ***
65  *
66  * Revision 2.0 1999/02/24 15:40:32 hyc
67  * *** empty log message ***
68  *
69  *
70  * $Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_math.C,v 1.6 2016/12/05 16:18:20 j_novak Exp $
71  *
72  */
73 
74 // Fichiers include
75 // ----------------
76 #include <cmath>
77 #include <cstdlib>
78 
79 #include "valeur.h"
80 
81 
82  //-------//
83  // Sinus //
84  //-------//
85 
86 namespace Lorene {
87 Valeur sin(const Valeur& ti)
88 {
89  // Protection
90  assert(ti.get_etat() != ETATNONDEF) ;
91 
92  // Cas ETATZERO
93  if (ti.get_etat() == ETATZERO) {
94  return ti ;
95  }
96 
97  // Cas general
98  assert(ti.get_etat() == ETATQCQ) ; // sinon...
99  if (ti.c == 0x0) { // Il faut la valeur physique
100  ti.coef_i() ;
101  }
102  Valeur to(ti.get_mg()) ; // Valeur resultat
103  to.set_etat_c_qcq() ;
104  *(to.c) = sin( *(ti.c) ) ;
105  return to ;
106 }
107 
108  //---------//
109  // Cosinus //
110  //---------//
111 
112 Valeur cos(const Valeur& ti)
113 {
114  // Protection
115  assert(ti.get_etat() != ETATNONDEF) ;
116 
117  Valeur to(ti.get_mg()) ; // Valeur resultat
118  to.set_etat_c_qcq() ;
119 
120  // Cas ETATZERO
121  if (ti.get_etat() == ETATZERO) {
122  *(to.c) = 1. ;
123  return to ;
124  }
125 
126  // Cas general
127  assert(ti.get_etat() == ETATQCQ) ; // sinon...
128  if (ti.c == 0x0) { // Il faut la valeur physique
129  ti.coef_i() ;
130  }
131  *(to.c) = cos( *(ti.c) ) ;
132  return to ;
133 }
134 
135  //----------//
136  // Tangente //
137  //----------//
138 
139 Valeur tan(const Valeur& ti)
140 {
141  // Protection
142  assert(ti.get_etat() != ETATNONDEF) ;
143 
144  // Cas ETATZERO
145  if (ti.get_etat() == ETATZERO) {
146  return ti ;
147  }
148 
149  // Cas general
150  assert(ti.get_etat() == ETATQCQ) ; // sinon...
151  if (ti.c == 0x0) { // Il faut la valeur physique
152  ti.coef_i() ;
153  }
154  Valeur to(ti.get_mg()) ; // Valeur resultat
155  to.set_etat_c_qcq() ;
156  *(to.c) = tan( *(ti.c) ) ;
157  return to ;
158 }
159 
160  //----------//
161  // ArcSinus //
162  //----------//
163 
164 Valeur asin(const Valeur& ti)
165 {
166  // Protection
167  assert(ti.get_etat() != ETATNONDEF) ;
168 
169  // Cas ETATZERO
170  if (ti.get_etat() == ETATZERO) {
171  return ti ;
172  }
173 
174  // Cas general
175  assert(ti.get_etat() == ETATQCQ) ; // sinon...
176  if (ti.c == 0x0) { // Il faut la valeur physique
177  ti.coef_i() ;
178  }
179  Valeur to(ti.get_mg()) ; // Valeur resultat
180  to.set_etat_c_qcq() ;
181  *(to.c) = asin( *(ti.c) ) ;
182  return to ;
183 }
184 
185  //------------//
186  // ArcCosinus //
187  //------------//
188 
189 Valeur acos(const Valeur& ti)
190 {
191  // Protection
192  assert(ti.get_etat() != ETATNONDEF) ;
193 
194  Valeur to(ti.get_mg()) ; // Valeur resultat
195  to.set_etat_c_qcq() ;
196 
197  // Cas ETATZERO
198  if (ti.get_etat() == ETATZERO) {
199  *(to.c) = M_PI * .5 ;
200  return to ;
201  }
202 
203  // Cas general
204  assert(ti.get_etat() == ETATQCQ) ; // sinon...
205  if (ti.c == 0x0) { // Il faut la valeur physique
206  ti.coef_i() ;
207  }
208  *(to.c) = acos( *(ti.c) ) ;
209  return to ;
210 }
211 
212  //-------------//
213  // ArcTangente //
214  //-------------//
215 
216 Valeur atan(const Valeur& ti)
217 {
218  // Protection
219  assert(ti.get_etat() != ETATNONDEF) ;
220 
221  // Cas ETATZERO
222  if (ti.get_etat() == ETATZERO) {
223  return ti ;
224  }
225 
226  // Cas general
227  assert(ti.get_etat() == ETATQCQ) ; // sinon...
228  if (ti.c == 0x0) { // Il faut la valeur physique
229  ti.coef_i() ;
230  }
231  Valeur to(ti.get_mg()) ; // Valeur resultat
232  to.set_etat_c_qcq() ;
233  *(to.c) = atan( *(ti.c) ) ;
234  return to ;
235 }
236 
237  //------//
238  // Sqrt //
239  //------//
240 
241 Valeur sqrt(const Valeur& ti)
242 {
243  // Protection
244  assert(ti.get_etat() != ETATNONDEF) ;
245 
246  // Cas ETATZERO
247  if (ti.get_etat() == ETATZERO) {
248  return ti ;
249  }
250 
251  // Cas general
252  assert(ti.get_etat() == ETATQCQ) ; // sinon...
253  if (ti.c == 0x0) { // Il faut la valeur physique
254  ti.coef_i() ;
255  }
256  Valeur to(ti.get_mg()) ; // Valeur resultat
257  to.set_etat_c_qcq() ;
258  *(to.c) = sqrt( *(ti.c) ) ;
259  return to ;
260 }
261 
262  //---------------//
263  // Exponantielle //
264  //---------------//
265 
266 Valeur exp(const Valeur& ti)
267 {
268  // Protection
269  assert(ti.get_etat() != ETATNONDEF) ;
270 
271  Valeur to(ti.get_mg()) ; // Valeur resultat
272  to.set_etat_c_qcq() ;
273 
274  // Cas ETATZERO
275  if (ti.get_etat() == ETATZERO) {
276  *(to.c) = 1. ;
277  return to ;
278  }
279 
280  // Cas general
281  assert(ti.get_etat() == ETATQCQ) ; // sinon...
282  if (ti.c == 0x0) { // Il faut la valeur physique
283  ti.coef_i() ;
284  }
285  *(to.c) = exp( *(ti.c) ) ;
286  return to ;
287 }
288 
289  //--------------------//
290  // Heaviside Function //
291  //--------------------//
292 
293 Valeur Heaviside(const Valeur& ti)
294 {
295  // Protection
296  assert(ti.get_etat() != ETATNONDEF) ;
297 
298  Valeur to(ti.get_mg()) ; // Valeur resultat
299  to.set_etat_c_qcq() ;
300 
301  // Cas ETATZERO
302  if (ti.get_etat() == ETATZERO) {
303  *(to.c) = 0. ;
304  return to ;
305  }
306 
307  // Cas general
308  assert(ti.get_etat() == ETATQCQ) ; // otherwise
309  if (ti.c == 0x0) { // Use the physical value << What?? (used Google translate)
310  ti.coef_i() ;
311  }
312 
313  *(to.c) = Heaviside( *(ti.c) ) ;
314 
315  return to ;
316 }
317 
318  //-------------//
319  // Log naturel //
320  //-------------//
321 
322 Valeur log(const Valeur& ti)
323 {
324  // Protection
325  assert(ti.get_etat() != ETATNONDEF) ;
326 
327  // Cas ETATZERO
328  if (ti.get_etat() == ETATZERO) {
329  cout << "Valeur log: log(ETATZERO) !" << endl ;
330  abort () ;
331  }
332 
333  // Cas general
334  assert(ti.get_etat() == ETATQCQ) ; // sinon...
335  if (ti.c == 0x0) { // Il faut la valeur physique
336  ti.coef_i() ;
337  }
338  Valeur to(ti.get_mg()) ; // Valeur resultat
339  to.set_etat_c_qcq() ;
340  *(to.c) = log( *(ti.c) ) ;
341  return to ;
342 }
343 
344  //-------------//
345  // Log decimal //
346  //-------------//
347 
348 Valeur log10(const Valeur& ti)
349 {
350  // Protection
351  assert(ti.get_etat() != ETATNONDEF) ;
352 
353  // Cas ETATZERO
354  if (ti.get_etat() == ETATZERO) {
355  cout << "Valeur log10: log10(ETATZERO) !" << endl ;
356  abort () ;
357  }
358 
359  // Cas general
360  assert(ti.get_etat() == ETATQCQ) ; // sinon...
361  if (ti.c == 0x0) { // Il faut la valeur physique
362  ti.coef_i() ;
363  }
364  Valeur to(ti.get_mg()) ; // Valeur resultat
365  to.set_etat_c_qcq() ;
366  *(to.c) = log10( *(ti.c) ) ;
367  return to ;
368 }
369 
370  //--------------//
371  // Power entier //
372  //--------------//
373 
374 Valeur pow(const Valeur& ti, int n)
375 {
376  // Protection
377  assert(ti.get_etat() != ETATNONDEF) ;
378 
379  // Cas ETATZERO
380  if (ti.get_etat() == ETATZERO) {
381  if (n > 0) {
382  return ti ;
383  }
384  else {
385  cout << "Valeur pow: ETATZERO^n with n<=0 ! "<< endl ;
386  abort () ;
387  }
388  }
389 
390  // Cas general
391  assert(ti.get_etat() == ETATQCQ) ; // sinon...
392  if (ti.c == 0x0) { // Il faut la valeur physique
393  ti.coef_i() ;
394  }
395  Valeur to(ti.get_mg()) ; // Valeur resultat
396  to.set_etat_c_qcq() ;
397  double x = n ;
398  *(to.c) = pow( *(ti.c), x ) ;
399  return to ;
400 }
401 
402  //--------------//
403  // Power double //
404  //--------------//
405 
406 Valeur pow(const Valeur& ti, double x)
407 {
408  // Protection
409  assert(ti.get_etat() != ETATNONDEF) ;
410 
411  // Cas ETATZERO
412  if (ti.get_etat() == ETATZERO) {
413  if (x > 0) {
414  return ti ;
415  }
416  else {
417  cout << "Valeur pow: ETATZERO^x with x<=0 !" << endl ;
418  abort () ;
419  }
420  }
421 
422  // Cas general
423  assert(ti.get_etat() == ETATQCQ) ; // sinon...
424  if (ti.c == 0x0) { // Il faut la valeur physique
425  ti.coef_i() ;
426  }
427  Valeur to(ti.get_mg()) ; // Valeur resultat
428  to.set_etat_c_qcq() ;
429  *(to.c) = pow( *(ti.c), x ) ;
430  return to ;
431 }
432 
433  //----------------//
434  // Absolute value //
435  //----------------//
436 
437 Valeur abs(const Valeur& vi)
438 {
439  // Protection
440  assert(vi.get_etat() != ETATNONDEF) ;
441 
442  // Cas ETATZERO
443  if (vi.get_etat() == ETATZERO) {
444  return vi ;
445  }
446 
447  // Cas general
448 
449  assert(vi.get_etat() == ETATQCQ) ; // sinon...
450  if (vi.c == 0x0) { // Il faut la valeur physique
451  vi.coef_i() ;
452  }
453 
454  Valeur vo(vi.get_mg()) ; // Valeur resultat
455 
456  vo.set_etat_c_qcq() ;
457 
458  *(vo.c) = abs( *(vi.c) ) ; // abs(Mtbl)
459 
460  return vo ;
461 }
462  //----------------//
463  // Cube root //
464  //----------------//
465 
466 Valeur racine_cubique(const Valeur& vi)
467 {
468 // Protection
469  assert(vi.get_etat() != ETATNONDEF) ;
470 
471  // Cas ETATZERO
472  if (vi.get_etat() == ETATZERO) {
473  return vi ;
474  }
475 
476  // Cas general
477  assert(vi.get_etat() == ETATQCQ) ; // sinon...
478  if (vi.c == 0x0) { // Il faut la valeur physique
479  vi.coef_i() ;
480  }
481  Valeur vo(vi.get_mg()) ; // Valeur resultat
482  vo.set_etat_c_qcq() ;
483  *(vo.c) = racine_cubique( *(vi.c) ) ;
484  return vo ;
485 }
486  //-------------------------------//
487  // totalmax //
488  //-------------------------------//
489 
490 double totalmax(const Valeur& vi) {
491 
492  // Protection
493  assert(vi.get_etat() != ETATNONDEF) ;
494 
495 // Tbl resu(vi.get_mg()->get_nzone()) ;
496  double resu ;
497 
498  if (vi.get_etat() == ETATZERO) {
499  resu = 0 ;
500  }
501  else {
502 
503  assert(vi.get_etat() == ETATQCQ) ;
504  if (vi.c == 0x0) { // Il faut la valeur physique
505  vi.coef_i() ;
506  }
507 
508  resu = totalmax( *(vi.c) ) ; // max(Mtbl)
509 
510  }
511 
512  return resu ;
513 }
514 
515  //-------------------------------//
516  // totalmin //
517  //-------------------------------//
518 
519 double totalmin(const Valeur& vi) {
520 
521  // Protection
522  assert(vi.get_etat() != ETATNONDEF) ;
523 
524  double resu ;
525 
526  if (vi.get_etat() == ETATZERO) {
527  resu = 0 ;
528  }
529  else {
530 
531  assert(vi.get_etat() == ETATQCQ) ;
532  if (vi.c == 0x0) { // Il faut la valeur physique
533  vi.coef_i() ;
534  }
535 
536  resu = totalmin( *(vi.c) ) ; // min(Mtbl)
537 
538  }
539 
540  return resu ;
541 }
542 
543  //-------------------------------//
544  // max //
545  //-------------------------------//
546 
547 Tbl max(const Valeur& vi) {
548 
549  // Protection
550  assert(vi.get_etat() != ETATNONDEF) ;
551 
552  Tbl resu(vi.get_mg()->get_nzone()) ;
553 
554  if (vi.get_etat() == ETATZERO) {
555  resu.annule_hard() ;
556  }
557  else {
558 
559  assert(vi.get_etat() == ETATQCQ) ;
560  if (vi.c == 0x0) { // Il faut la valeur physique
561  vi.coef_i() ;
562  }
563 
564  resu = max( *(vi.c) ) ; // max(Mtbl)
565 
566  }
567 
568  return resu ;
569 }
570 
571  //-------------------------------//
572  // min //
573  //-------------------------------//
574 
575 Tbl min(const Valeur& vi) {
576 
577  // Protection
578  assert(vi.get_etat() != ETATNONDEF) ;
579 
580  Tbl resu(vi.get_mg()->get_nzone()) ;
581 
582  if (vi.get_etat() == ETATZERO) {
583  resu.annule_hard() ;
584  }
585  else {
586 
587  assert(vi.get_etat() == ETATQCQ) ;
588  if (vi.c == 0x0) { // Il faut la valeur physique
589  vi.coef_i() ;
590  }
591 
592  resu = min( *(vi.c) ) ; // min(Mtbl)
593 
594  }
595 
596  return resu ;
597 }
598 
599  //-------------------------------//
600  // norme //
601  //-------------------------------//
602 
603 Tbl norme(const Valeur& vi) {
604 
605  // Protection
606  assert(vi.get_etat() != ETATNONDEF) ;
607 
608  Tbl resu(vi.get_mg()->get_nzone()) ;
609 
610  if (vi.get_etat() == ETATZERO) {
611  resu.annule_hard() ;
612  }
613  else {
614 
615  assert(vi.get_etat() == ETATQCQ) ;
616  if (vi.c == 0x0) { // Il faut la valeur physique
617  vi.coef_i() ;
618  }
619 
620  resu = norme( *(vi.c) ) ; // norme(Mtbl)
621 
622  }
623 
624  return resu ;
625 }
626 
627  //-------------------------------//
628  // diffrel //
629  //-------------------------------//
630 
631 Tbl diffrel(const Valeur& v1, const Valeur& v2) {
632 
633  // Protections
634  assert(v1.get_etat() != ETATNONDEF) ;
635  assert(v2.get_etat() != ETATNONDEF) ;
636 
637  int nz = v1.get_mg()->get_nzone() ;
638  Tbl resu(nz) ;
639 
640  Valeur diff = v1 - v2 ;
641  Tbl normdiff = norme(diff) ;
642  Tbl norme2 = norme(v2) ;
643 
644  assert(normdiff.get_etat() == ETATQCQ) ;
645  assert(norme2.get_etat() == ETATQCQ) ;
646 
647  resu.set_etat_qcq() ;
648  for (int l=0; l<nz; l++) {
649  if ( norme2(l) == double(0) ) {
650  resu.set(l) = normdiff(l) ;
651  }
652  else{
653  resu.set(l) = normdiff(l) / norme2(l) ;
654  }
655  }
656 
657  return resu ;
658 
659 }
660 
661  //-------------------------------//
662  // diffrelmax //
663  //-------------------------------//
664 
665 Tbl diffrelmax(const Valeur& v1, const Valeur& v2) {
666 
667  // Protections
668  assert(v1.get_etat() != ETATNONDEF) ;
669  assert(v2.get_etat() != ETATNONDEF) ;
670 
671  int nz = v1.get_mg()->get_nzone() ;
672  Tbl resu(nz) ;
673 
674  Tbl max2 = max(abs(v2)) ;
675  Valeur diff = v1 - v2 ;
676  Tbl maxdiff = max(abs(diff)) ;
677 
678  assert(maxdiff.get_etat() == ETATQCQ) ;
679  assert(max2.get_etat() == ETATQCQ) ;
680 
681  resu.set_etat_qcq() ;
682  for (int l=0; l<nz; l++) {
683  if ( max2(l) == double(0) ) {
684  resu.set(l) = maxdiff(l) ;
685  }
686  else{
687  resu.set(l) = maxdiff(l) / max2(l) ;
688  }
689  }
690 
691  return resu ;
692 
693 }
694 
695 }
Lorene::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
Lorene::asin
Cmp asin(const Cmp &)
Arcsine.
Definition: cmp_math.C:147
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:273
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
Lorene
Lorene prototypes.
Definition: app_hor.h:67
Lorene::racine_cubique
Cmp racine_cubique(const Cmp &)
Cube root.
Definition: cmp_math.C:248
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Lorene::Valeur::coef_i
void coef_i() const
Computes the physical value of *this.
Definition: valeur_coef_i.C:143
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
Lorene::Tbl::get_etat
int get_etat() const
Gives the logical state.
Definition: tbl.h:414
Lorene::min
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition: cmp_math.C:461
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:97
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:507
Lorene::Valeur::get_mg
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
Definition: valeur.h:763
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
Lorene::tan
Cmp tan(const Cmp &)
Tangent.
Definition: cmp_math.C:123
Lorene::Valeur::get_etat
int get_etat() const
Returns the logical state.
Definition: valeur.h:760
Lorene::norme
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:484
Lorene::atan
Cmp atan(const Cmp &)
Arctangent.
Definition: cmp_math.C:198
Lorene::Valeur::c
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:309
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:438
Lorene::totalmax
double totalmax(const Mtbl &)
Maximum value of the Mtbl elements in all domains.
Definition: mtbl_math.C:497
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Lorene::totalmin
double totalmin(const Mtbl &)
Minimum value of the Mtbl elements in all domain.
Definition: mtbl_math.C:525
Lorene::log10
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:325
Lorene::acos
Cmp acos(const Cmp &)
Arccosine.
Definition: cmp_math.C:172
Lorene::abs
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:413
Lorene::Valeur::set_etat_c_qcq
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
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:72
Lorene::Heaviside
Mtbl Heaviside(const Mtbl &)
Heaviside function.
Definition: mtbl_math.C:320
Lorene::Tbl
Basic array class.
Definition: tbl.h:164
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:375
Lorene::diffrelmax
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:542