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
86namespace Lorene {
87Valeur 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
112Valeur 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
139Valeur 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
164Valeur 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
189Valeur 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
216Valeur 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
241Valeur 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
266Valeur 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
293Valeur 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
322Valeur 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
348Valeur 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
374Valeur 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
406Valeur 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
437Valeur 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
466Valeur 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
490double 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
519double 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
547Tbl 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
575Tbl 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
603Tbl 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
631Tbl 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
665Tbl 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::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:465
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::get_etat
int get_etat() const
Gives the logical state.
Definition tbl.h:394
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:375
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:364
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:297
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::Valeur::get_etat
int get_etat() const
Returns the logical state.
Definition valeur.h:760
Lorene::Valeur::c
Mtbl * c
Values of the function at the points of the multi-grid.
Definition valeur.h:309
Lorene::Valeur::coef_i
void coef_i() const
Computes the physical value of *this.
Definition valeur_coef_i.C:143
Lorene::Valeur::get_mg
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
Definition valeur.h:763
Lorene::atan
Cmp atan(const Cmp &)
Arctangent.
Definition cmp_math.C:198
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:223
Lorene::log10
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition cmp_math.C:325
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:273
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:72
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:507
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::acos
Cmp acos(const Cmp &)
Arccosine.
Definition cmp_math.C:172
Lorene::asin
Cmp asin(const Cmp &)
Arcsine.
Definition cmp_math.C:147
Lorene::racine_cubique
Cmp racine_cubique(const Cmp &)
Cube root.
Definition cmp_math.C:248
Lorene::min
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition cmp_math.C:461
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:438
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:351
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:97
Lorene::abs
Cmp abs(const Cmp &)
Absolute value.
Definition cmp_math.C:413
Lorene::diffrelmax
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition cmp_math.C:542
Lorene::tan
Cmp tan(const Cmp &)
Tangent.
Definition cmp_math.C:123
Lorene::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:299
Lorene::Heaviside
Mtbl Heaviside(const Mtbl &)
Heaviside function.
Definition mtbl_math.C:320
Lorene::totalmin
double totalmin(const Mtbl &)
Minimum value of the Mtbl elements in all domain.
Definition mtbl_math.C:525
Lorene::totalmax
double totalmax(const Mtbl &)
Maximum value of the Mtbl elements in all domains.
Definition mtbl_math.C:497
Lorene
Lorene prototypes.
Definition app_hor.h:67
Lorene::x
Coord x
x coordinate centered on the grid
Definition map.h:738