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