LORENE
valeur_arithm.C
1 /*
2  * Arithmetical operations for class Valeur
3  *
4  */
5 
6 /*
7  * Copyright (c) 1999-2000 Jean-Alain Marck
8  * Copyright (c) 1999-2005 Eric Gourgoulhon
9  * Copyright (c) 1999-2001 Philippe Grandclement
10  * Copyright (c) 2004 Jerome Novak
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License as published by
16  * the Free Software Foundation; either version 2 of the License, or
17  * (at your option) any later version.
18  *
19  * LORENE is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with LORENE; if not, write to the Free Software
26  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27  *
28  */
29 
30 
31 char valeur_arithm_C[] = "$Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_arithm.C,v 1.7 2014/10/13 08:53:49 j_novak Exp $" ;
32 
33 /*
34  * $Id: valeur_arithm.C,v 1.7 2014/10/13 08:53:49 j_novak Exp $
35  * $Log: valeur_arithm.C,v $
36  * Revision 1.7 2014/10/13 08:53:49 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.6 2014/10/06 15:13:22 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.5 2008/08/27 08:52:55 jl_cornou
43  * Added Jacobi(0,2) polynomials case
44  *
45  * Revision 1.4 2005/11/17 15:19:23 e_gourgoulhon
46  * Added Valeur + Mtbl and Valeur - Mtbl.
47  *
48  * Revision 1.3 2004/07/06 13:36:30 j_novak
49  * Added methods for desaliased product (operator |) only in r direction.
50  *
51  * Revision 1.2 2002/10/16 14:37:15 j_novak
52  * Reorganization of #include instructions of standard C++, in order to
53  * use experimental version 3 of gcc.
54  *
55  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
56  * LORENE
57  *
58  * Revision 2.20 2001/08/31 15:23:48 eri * operator% : traitement du cas ou le Tbl est zero dans une zone.
59  *
60  * Revision 2.19 2001/05/28 12:42:27 eric
61  * Passage en ylm pour le desaliasing dans operator%.
62  *
63  * Revision 2.18 2001/05/26 14:50:21 eric
64  * Ajout de l'operator% : produit de deux Valeur avec desaliasage.
65  *
66  * Revision 2.17 2000/01/14 14:42:47 eric
67  * Valeur * double : operation effectue dans l'espace des coefficients si
68  * la Valeur n'est pas a jour ds l'espace des config.
69  * Valeur / double : idem.
70  * += Valeur : idem.
71  * -= Valeur : idem.
72  *
73  * Pour += Valeur et -= Valeur les schemas sont desormais calques sur
74  * Valeur + Valeur et Valeur - Valeur.
75  *
76  * Revision 2.16 1999/12/10 16:33:36 eric
77  * Dans l'arithmetique membre (+=, -=, *=), on n'appelle desormais
78  * del_deriv() que tout a la fin.
79  *
80  * Revision 2.15 1999/11/30 14:12:54 phil
81  * gestion de base dans operator/ (double,Vlaeur)
82  *
83  * Revision 2.14 1999/11/30 12:42:10 eric
84  * Le membre base est desormais un objet de type Base_val et non plus
85  * un pointeur vers une Base_val.
86  *
87  * Revision 2.13 1999/11/29 13:28:06 eric
88  * *** empty log message ***
89  *
90  * Revision 2.12 1999/11/29 10:20:50 eric
91  * Ajout de Valeur/Mtbl et Mtbl / Valeur.
92  *
93  * Revision 2.11 1999/11/29 10:06:05 eric
94  * Ajout de Valeur*Mtbl et Mtbl*Valeur
95  *
96  * Revision 2.10 1999/10/26 14:40:29 phil
97  * On gere les bases pour *, *=, /= et /
98  *
99  * Revision 2.9 1999/10/20 15:31:28 eric
100  * Ajout de la plupart des fonctions arithmetiques.
101  *
102  * Revision 2.8 1999/10/18 15:07:47 eric
103  * La fonction membre annule() est rebaptisee annule_hard().
104  *
105  * Revision 2.7 1999/10/13 15:50:56 eric
106  * Depoussierage.
107  *
108  * Revision 2.6 1999/09/15 10:02:26 phil
109  * gestion de la base dans Valeur operator (double, const Valeur &)
110  *
111  * Revision 2.5 1999/09/14 17:18:36 phil
112  * aout de Valeur operator* (double, const Valeur&)
113  *
114  * Revision 2.4 1999/09/13 14:52:55 phil
115  * *** empty log message ***
116  *
117  * Revision 2.3 1999/04/09 12:38:05 phil
118  * *** empty log message ***
119  *
120  * Revision 2.2 1999/04/09 12:26:59 phil
121  * Ajout de valeur * coord
122  *
123  * Revision 2.1 1999/02/22 15:49:28 hyc
124  * *** empty log message ***
125  *
126  *
127  * $Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_arithm.C,v 1.7 2014/10/13 08:53:49 j_novak Exp $
128  *
129  */
130 
131 // Fichiers include
132 // ----------------
133 #include <cmath>
134 #include <cassert>
135 #include <cstdlib>
136 
137 #include "mtbl.h"
138 #include "valeur.h"
139 #include "coord.h"
140  //********************//
141  // OPERATEURS UNAIRES //
142  //********************//
143 
144 namespace Lorene {
145 Valeur operator+(const Valeur & vi) {
146 
147  // Protection
148  assert(vi.get_etat() != ETATNONDEF) ;
149 
150  return vi ;
151 }
152 
153 Valeur operator-(const Valeur & vi) {
154 
155  // Protection
156  assert(vi.get_etat() != ETATNONDEF) ;
157 
158  // Cas particulier
159  if (vi.get_etat() == ETATZERO) {
160  return vi ;
161  }
162 
163  // Cas general
164  assert(vi.get_etat() == ETATQCQ) ; // sinon...
165 
166  Valeur resu(vi.get_mg()) ; // Valeur resultat
167 
168  if (vi.c != 0x0) {
169  resu = - *(vi.c) ;
170  resu.base = vi.base ; // N'oublions pas la base...
171  }
172  else{
173  assert(vi.c_cf != 0x0) ;
174  resu = - *(vi.c_cf) ; // Dans ce cas la base est prise en
175  // charge par l'operator=(const Mtbl_cf&).
176  }
177 
178  // Termine
179  return resu ;
180 }
181 
182  //**********//
183  // ADDITION //
184  //**********//
185 
186 // Valeur + Valeur
187 // ---------------
188 Valeur operator+(const Valeur& t1, const Valeur& t2)
189 {
190  // Protections
191  assert(t1.get_etat() != ETATNONDEF) ;
192  assert(t2.get_etat() != ETATNONDEF) ;
193  assert(t1.get_mg() == t2.get_mg()) ;
194 
195  // Cas particuliers
196  if (t1.get_etat() == ETATZERO) {
197  return t2 ;
198  }
199  if (t2.get_etat() == ETATZERO) {
200  return t1 ;
201  }
202  assert(t1.get_etat() == ETATQCQ) ; // sinon...
203  assert(t2.get_etat() == ETATQCQ) ; // sinon...
204 
205  Valeur resu(t1.get_mg()) ;
206 
207  // On privelegie l'addition dans l'espace des configurations:
208  // ----------------------------------------------------------
209  if (t1.c != 0x0) {
210  if (t2.c != 0x0) {
211  resu = *(t1.c) + *(t2.c) ; // Addition des Mtbl
212  resu.base = t1.base ;
213  }
214  else {
215  assert(t2.c_cf != 0x0) ;
216  if (t1.c_cf != 0x0) {
217  resu = *(t1.c_cf) + *(t2.c_cf) ; // Addition des Mtbl_cf
218  }
219  else {
220  t2.coef_i() ;
221  resu = *(t1.c) + *(t2.c) ; // Addition des Mtbl
222  resu.base = t1.base ;
223  }
224  }
225  }
226  else{ // Cas ou t1.c n'est pas a jour
227  assert(t1.c_cf != 0x0) ;
228  if (t2.c_cf != 0x0) {
229  resu = *(t1.c_cf) + *(t2.c_cf) ; // Addition des Mtbl_cf
230  }
231  else {
232  assert(t2.c != 0x0) ;
233  t1.coef_i() ;
234  resu = *(t1.c) + *(t2.c) ; // Addition des Mtbl
235  resu.base = t1.base ;
236  }
237  }
238 
239  return resu ;
240 }
241 
242 // Valeur + Mtbl
243 // -------------
244 Valeur operator+(const Valeur& t1, const Mtbl& mi)
245 {
246  // Protections
247  assert(t1.get_etat() != ETATNONDEF) ;
248  assert(mi.get_etat() != ETATNONDEF) ;
249 
250  // Cas particuliers
251  if (mi.get_etat() == ETATZERO) {
252  return t1 ;
253  }
254 
255  Valeur resu(t1) ;
256 
257  if (t1.get_etat() == ETATZERO) {
258  resu.set_etat_c_qcq() ;
259  *(resu.c) = mi ;
260  }
261  else{
262  assert(resu.get_etat() == ETATQCQ) ; // sinon ...
263  resu.coef_i() ; // l'addition se fait dans l'espace des configurations
264  resu.set_etat_c_qcq() ;
265  *(resu.c) += mi ;
266  }
267 
268  return resu ;
269 }
270 
271 // Mtbl + Valeur
272 // -------------
273 Valeur operator+(const Mtbl& mi, const Valeur& t1) {
274  return t1 + mi ;
275 }
276 
277 
278 // Valeur + double
279 // ---------------
280 Valeur operator+(const Valeur& t1, double x)
281 {
282  // Protections
283  assert(t1.get_etat() != ETATNONDEF) ;
284 
285  // Cas particuliers
286  if (x == double(0)) {
287  return t1 ;
288  }
289 
290  Valeur resu(t1) ;
291 
292  if (t1.get_etat() == ETATZERO) {
293  resu.set_etat_c_qcq() ;
294  *(resu.c) = x ;
295  }
296  else{
297  assert(resu.get_etat() == ETATQCQ) ; // sinon ...
298  resu.coef_i() ; // l'addition se fait dans l'espace des configurations
299  resu.set_etat_c_qcq() ;
300  *(resu.c) = *(resu.c) + x ;
301  }
302 
303  return resu ;
304 }
305 
306 // double + Valeur
307 // ---------------
308 Valeur operator+(double x, const Valeur& t1) // double + Valeur
309 {
310  return t1 + x ;
311 }
312 
313 // Valeur + int
314 // ------------
315 Valeur operator+(const Valeur& t1, int m) // Valeur + int
316 {
317  return t1 + double(m) ;
318 }
319 
320 // int + Valeur
321 // -------------
322 Valeur operator+(int m, const Valeur& t1) // int + Valeur
323 {
324  return t1 + double(m) ;
325 }
326 
327 
328 
329  //**************//
330  // SOUSTRACTION //
331  //**************//
332 
333 // Valeur - Valeur
334 // ---------------
335 Valeur operator-(const Valeur& t1, const Valeur& t2)
336 {
337  // Protections
338  assert(t1.get_etat() != ETATNONDEF) ;
339  assert(t2.get_etat() != ETATNONDEF) ;
340  assert(t1.get_mg() == t2.get_mg()) ;
341 
342  // Cas particuliers
343  if (t1.get_etat() == ETATZERO) {
344  return - t2 ;
345  }
346  if (t2.get_etat() == ETATZERO) {
347  return t1 ;
348  }
349  assert(t1.get_etat() == ETATQCQ) ; // sinon...
350  assert(t2.get_etat() == ETATQCQ) ; // sinon...
351 
352  Valeur resu(t1.get_mg()) ;
353 
354  // On privelegie la soustraction dans l'espace des configurations:
355  // ---------------------------------------------------------------
356  if (t1.c != 0x0) {
357  if (t2.c != 0x0) {
358  resu = *(t1.c) - *(t2.c) ; // Soustraction des Mtbl
359  resu.base = t1.base ;
360  }
361  else {
362  assert(t2.c_cf != 0x0) ;
363  if (t1.c_cf != 0x0) {
364  resu = *(t1.c_cf) - *(t2.c_cf) ; // Soustraction des Mtbl_cf
365  }
366  else {
367  t2.coef_i() ;
368  resu = *(t1.c) - *(t2.c) ; // Soustraction des Mtbl
369  resu.base = t1.base ;
370  }
371  }
372  }
373  else{ // Cas ou t1.c n'est pas a jour
374  assert(t1.c_cf != 0x0) ;
375  if (t2.c_cf != 0x0) {
376  resu = *(t1.c_cf) - *(t2.c_cf) ; // Soustraction des Mtbl_cf
377  }
378  else {
379  assert(t2.c != 0x0) ;
380  t1.coef_i() ;
381  resu = *(t1.c) - *(t2.c) ; // Soustraction des Mtbl
382  resu.base = t1.base ;
383  }
384  }
385 
386  return resu ;
387 }
388 
389 
390 // Valeur - Mtbl
391 // -------------
392 Valeur operator-(const Valeur& t1, const Mtbl& mi)
393 {
394  // Protections
395  assert(t1.get_etat() != ETATNONDEF) ;
396  assert(mi.get_etat() != ETATNONDEF) ;
397 
398  // Cas particuliers
399  if (mi.get_etat() == ETATZERO) {
400  return t1 ;
401  }
402 
403  Valeur resu(t1) ;
404 
405  if (t1.get_etat() == ETATZERO) {
406  resu.set_etat_c_qcq() ;
407  *(resu.c) = - mi ;
408  }
409  else{
410  assert(resu.get_etat() == ETATQCQ) ; // sinon ...
411  resu.coef_i() ; // substraction in configuration space
412  resu.set_etat_c_qcq() ;
413  *(resu.c) -= mi ;
414  }
415 
416  return resu ;
417 }
418 
419 // Mtbl - Valeur
420 // -------------
421 Valeur operator-(const Mtbl& mi, const Valeur& t1) {
422  return - (t1 - mi) ;
423 }
424 
425 // Valeur - double
426 // ---------------
427 Valeur operator-(const Valeur& t1, double x)
428 {
429  // Protections
430  assert(t1.get_etat() != ETATNONDEF) ;
431 
432  // Cas particuliers
433  if (x == double(0)) {
434  return t1 ;
435  }
436 
437  Valeur resu(t1) ;
438 
439  if (t1.get_etat() == ETATZERO) {
440  resu.set_etat_c_qcq() ;
441  *(resu.c) = - x ;
442  }
443  else{
444  assert(resu.get_etat() == ETATQCQ) ; // sinon ...
445  resu.coef_i() ; // l'addition se fait dans l'espace des configurations
446  resu.set_etat_c_qcq() ;
447  *(resu.c) = *(resu.c) - x ;
448  }
449 
450  return resu ;
451 }
452 
453 // double - Valeur
454 // ---------------
455 Valeur operator-(double x, const Valeur& t1) // double - Valeur
456 {
457  return - (t1 - x) ;
458 }
459 
460 // Valeur - int
461 // ------------
462 Valeur operator-(const Valeur& t1, int m) // Valeur - int
463 {
464  return t1 - double(m) ;
465 }
466 
467 // int - Valeur
468 // -------------
469 Valeur operator-(int m, const Valeur& t1) // int - Valeur
470 {
471  return double(m) - t1 ;
472 }
473 
474 
475 
476 
477 
478  //****************//
479  // MULTIPLICATION //
480  //****************//
481 
482 // Valeur * Valeur
483 // ---------------
484 Valeur operator*(const Valeur& t1, const Valeur& t2)
485 {
486  // Protections
487  assert(t1.get_etat() != ETATNONDEF) ;
488  assert(t2.get_etat() != ETATNONDEF) ;
489  assert(t1.get_mg() == t2.get_mg()) ;
490 
491  // Cas particuliers
492  if (t1.get_etat() == ETATZERO) {
493  return t1 ;
494  }
495  if (t2.get_etat() == ETATZERO) {
496  return t2 ;
497  }
498  assert(t1.get_etat() == ETATQCQ) ; // sinon...
499  assert(t2.get_etat() == ETATQCQ) ; // sinon...
500 
501  Valeur resu(t1.get_mg()) ;
502 
503  // La multiplication est faite dans l'espace des configurations:
504  if (t1.c == 0x0) {
505  t1.coef_i() ;
506  }
507  if (t2.c == 0x0) {
508  t2.coef_i() ;
509  }
510 
511  resu = (*(t1.c)) * (*(t2.c)) ; // Multiplication des Mtbl
512 
513  // affectation de la base :
514  resu.base = t1.base * t2.base;
515 
516  return resu ;
517 }
518 
519 
520 // Valeur * double
521 // ---------------
522 Valeur operator*(const Valeur& c1, double a) {
523 
524  // Protection
525  assert(c1.get_etat() != ETATNONDEF) ;
526 
527  // Cas particulier
528  if ((c1.get_etat() == ETATZERO) || ( a == double(1) )) {
529  return c1 ;
530  }
531 
532  // Cas general
533  assert(c1.get_etat() == ETATQCQ) ; // sinon...
534 
535  Valeur result(c1.get_mg()) ;
536 
537  if (c1.c != 0x0) {
538  result = *(c1.c) * a ; // Mtbl * double
539  result.base = c1.base ; // in this case, result.base must be set
540  // by hand
541  }
542  else {
543  assert(c1.c_cf != 0x0) ;
544  result = *(c1.c_cf) * a ; // Mtbl_cf * double
545  }
546 
547  return result ;
548 }
549 
550 // double * Valeur
551 // ---------------
552 Valeur operator*(double x, const Valeur& t1) // double * Valeur
553 {
554  return t1 * x ;
555 }
556 
557 // Valeur * int
558 // ------------
559 Valeur operator*(const Valeur& t1, int m) // Valeur * int
560 {
561  return t1 * double(m) ;
562 }
563 
564 // int * Valeur
565 // ------------
566 Valeur operator*(int m, const Valeur& t1) // int * Valeur
567 {
568  return t1 * double(m) ;
569 }
570 
571 
572 // Valeur * Mtbl
573 // --------------
574 Valeur operator*(const Valeur & c1, const Mtbl& c2) {
575 
576  // Protection
577  assert(c1.get_etat() != ETATNONDEF) ;
578 
579  // Cas particulier
580  if (c1.get_etat() == ETATZERO) {
581  return c1 ;
582  }
583 
584  // Cas general
585 
586  assert(c1.get_etat() == ETATQCQ) ;
587 
588  Valeur result(c1.get_mg()) ;
589 
590  // La multiplication est faite dans l'espace des configurations:
591  if (c1.c == 0x0) {
592  c1.coef_i() ;
593  }
594  result = *(c1.c) * c2 ; // Mtbl * Mtbl
595 
596  return result ;
597 }
598 
599 // Mtbl * Valeur
600 // --------------
601 Valeur operator*(const Mtbl& c1, const Valeur& t1) {
602 
603  return t1 * c1 ;
604 
605 }
606 
607 
608 // Valeur * Coord
609 // --------------
610 Valeur operator*(const Valeur & c1, const Coord& c2) {
611 
612  // Protection
613  assert(c1.get_etat() != ETATNONDEF) ;
614 
615  // Cas particulier
616  if (c1.get_etat() == ETATZERO) {
617  return c1 ;
618  }
619 
620  // Cas general
621 
622  assert(c1.get_etat() == ETATQCQ) ;
623 
624  Valeur result(c1.get_mg()) ;
625 
626  // La multiplication est faite dans l'espace des configurations:
627  if (c1.c == 0x0) {
628  c1.coef_i() ;
629  }
630  result = *(c1.c) * c2 ; // Mtbl * Coord
631 
632  return result ;
633 }
634 
635 // Coord * Valeur
636 // --------------
637 Valeur operator*(const Coord& c1, const Valeur& t1) {
638 
639  return t1 * c1 ;
640 
641 }
642 
643  //**********//
644  // DIVISION //
645  //**********//
646 
647 // Valeur / Valeur
648 // ---------------
649 Valeur operator/(const Valeur& t1, const Valeur& t2)
650 {
651  // Protections
652  assert(t1.get_etat() != ETATNONDEF) ;
653  assert(t2.get_etat() != ETATNONDEF) ;
654  assert(t1.get_mg() == t2.get_mg()) ;
655 
656  // Cas particuliers
657  if (t2.get_etat() == ETATZERO) {
658  cout << "Division by 0 in Valeur / Valeur !" << endl ;
659  abort() ;
660  }
661  if (t1.get_etat() == ETATZERO) {
662  return t1 ;
663  }
664 
665  assert(t1.get_etat() == ETATQCQ) ; // sinon...
666  assert(t2.get_etat() == ETATQCQ) ; // sinon...
667 
668  Valeur resu(t1.get_mg()) ;
669 
670  // La division est faite dans l'espace des configurations:
671  if (t1.c == 0x0) {
672  t1.coef_i() ;
673  }
674  if (t2.c == 0x0) {
675  t2.coef_i() ;
676  }
677 
678  resu = (*(t1.c)) / (*(t2.c)) ; // Division des Mtbl
679 
680  // affectation de la base :
681  resu.base = t1.base * t2.base;
682 
683  return resu ;
684 }
685 
686 // Valeur / double
687 // ---------------
688 Valeur operator/(const Valeur& t1, double x)
689 {
690  // Protections
691  assert(t1.get_etat() != ETATNONDEF) ;
692 
693  // Cas particuliers
694  if ( x == double(0) ) {
695  cout << "Division by 0 in Valeur / double !" << endl ;
696  abort() ;
697  }
698  if ((t1.get_etat() == ETATZERO) || ( x == double(1) )) {
699  return t1 ;
700  }
701 
702  assert(t1.get_etat() == ETATQCQ) ; // sinon...
703 
704  Valeur resu(t1.get_mg()) ;
705 
706  if (t1.c != 0x0) {
707  resu = *(t1.c) / x ; // Mtbl / double
708  resu.base = t1.base ; // in this case, resu.base must be set by hand
709  }
710  else {
711  assert(t1.c_cf != 0x0) ;
712  resu = *(t1.c_cf) / x ; // Mtbl_cf * double
713  }
714 
715  return resu ;
716 }
717 
718 // double / Valeur
719 // ---------------
720 Valeur operator/(double x, const Valeur & c2) {
721 
722  // Protection
723  assert(c2.get_etat() != ETATNONDEF) ;
724 
725  // Cas particuliers
726  if (c2.get_etat() == ETATZERO) {
727  cout << "Division by 0 in double / Valeur !" << endl ;
728  abort() ;
729  }
730 
731  // Cas general
732  assert(c2.get_etat() == ETATQCQ) ; // sinon...
733 
734  // Il faut les valeurs physiques de c2
735  if (c2.c == 0x0) {
736  c2.coef_i() ;
737  }
738 
739  // Le resultat
740  Valeur r(c2.get_mg()) ;
741  r.set_etat_c_qcq() ;
742  *(r.c) = x / *(c2.c) ;
743 
744  // affectation de la base :
745  r.base = c2.get_mg()->std_base_scal() * c2.base ;
746 
747  // Termine
748  return r ;
749 }
750 
751 // Valeur / int
752 // ------------
753 Valeur operator/(const Valeur& t1, int m) {
754  return t1 / double(m) ;
755 }
756 
757 // int / Valeur
758 // ------------
759 Valeur operator/(int m, const Valeur& t1) {
760  return double(m) / t1 ;
761 }
762 
763 // Valeur / Mtbl
764 // ---------------
765 Valeur operator/(const Valeur& t1, const Mtbl& m2)
766 {
767  // Protections
768  assert(t1.get_etat() != ETATNONDEF) ;
769  assert(m2.get_etat() != ETATNONDEF) ;
770 
771  // Cas particuliers
772  if ( m2.get_etat() == ETATZERO ) {
773  cout << "Division by 0 in Valeur / Mtbl !" << endl ;
774  abort() ;
775  }
776  if (t1.get_etat() == ETATZERO) {
777  return t1 ;
778  }
779 
780  assert(t1.get_etat() == ETATQCQ) ; // sinon...
781 
782  Valeur resu(t1.get_mg()) ;
783 
784  // La division est faite dans l'espace des configurations:
785  if (t1.c == 0x0) {
786  t1.coef_i() ;
787  }
788 
789  resu = (*(t1.c)) / m2 ; // Division Mtbl / Mtbl
790 
791  return resu ;
792 }
793 
794 // Mtbl / Valeur
795 // ---------------
796 Valeur operator/(const Mtbl& m1, const Valeur & c2) {
797 
798  // Protection
799  assert(m1.get_etat() != ETATNONDEF) ;
800  assert(c2.get_etat() != ETATNONDEF) ;
801 
802  // Cas particuliers
803  if (c2.get_etat() == ETATZERO) {
804  cout << "Division by 0 in Mtbl / Valeur !" << endl ;
805  abort() ;
806  }
807 
808  // Cas general
809  assert(c2.get_etat() == ETATQCQ) ; // sinon...
810 
811  // Le resultat
812  Valeur resu(c2.get_mg()) ;
813 
814  // Il faut les valeurs physiques de c2
815  if (c2.c == 0x0) {
816  c2.coef_i() ;
817  }
818 
819  resu = m1 / (*(c2.c)) ; // Division Mtbl / Mtbl
820 
821  // Termine
822  return resu ;
823 }
824 
825 
826 
827 
828 
829 
830  //*******************//
831  // operateurs +=,... //
832  //*******************//
833 
834 void Valeur::operator+=(const Valeur & vi) {
835 
836  // Protection
837  assert(mg == vi.get_mg()) ; // meme grille
838  assert(etat != ETATNONDEF) ; // etat defini
839  assert(vi.get_etat() != ETATNONDEF) ; // etat defini
840 
841  // Cas particulier
842  if (vi.get_etat() == ETATZERO) {
843  return ;
844  }
845 
846  // Cas general
847 
848  // Cas de l'etat ZERO
849  if (etat == ETATZERO) {
850  annule_hard() ;
851  }
852 
853 
854  if (c != 0x0) {
855  if (vi.c != 0x0) {
856  *c += *(vi.c) ; // += Mtbl
857  delete c_cf ;
858  c_cf = 0x0 ;
859  }
860  else {
861  assert(vi.c_cf != 0x0) ;
862  if (c_cf != 0x0) {
863  *c_cf += *(vi.c_cf) ; // += Mtbl_cf
864  delete c ;
865  c = 0x0 ;
866  }
867  else {
868  vi.coef_i() ;
869  *c += *(vi.c) ; // += Mtbl
870  delete c_cf ;
871  c_cf = 0x0 ;
872  }
873  }
874  }
875  else{ // Case where c is not up to date
876  assert(c_cf != 0x0) ;
877  if (vi.c_cf != 0x0) {
878  *c_cf += *(vi.c_cf) ; // += Mtbl_cf
879  delete c ;
880  c = 0x0 ;
881  }
882  else {
883  assert(vi.c != 0x0) ;
884  coef_i() ;
885  *c += *(vi.c) ; // += Mtbl
886  delete c_cf ;
887  c_cf = 0x0 ;
888  }
889  }
890 
891  // Menage (a ne faire qu'a la fin seulement)
892  del_deriv() ;
893 
894 
895 }
896 
897 void Valeur::operator-=(const Valeur & vi) {
898 
899  // Protection
900  assert(mg == vi.get_mg()) ; // meme grille
901  assert(etat != ETATNONDEF) ; // etat defini
902  assert(vi.get_etat() != ETATNONDEF) ; // etat defini
903 
904  // Cas particulier
905  if (vi.get_etat() == ETATZERO) {
906  return ;
907  }
908 
909  // Cas general
910 
911  // Cas de l'etat ZERO
912  if (etat == ETATZERO) {
913  annule_hard() ;
914  }
915 
916  if (c != 0x0) {
917  if (vi.c != 0x0) {
918  *c -= *(vi.c) ; // -= Mtbl
919  delete c_cf ;
920  c_cf = 0x0 ;
921  }
922  else {
923  assert(vi.c_cf != 0x0) ;
924  if (c_cf != 0x0) {
925  *c_cf -= *(vi.c_cf) ; // -= Mtbl_cf
926  delete c ;
927  c = 0x0 ;
928  }
929  else {
930  vi.coef_i() ;
931  *c -= *(vi.c) ; // -= Mtbl
932  delete c_cf ;
933  c_cf = 0x0 ;
934  }
935  }
936  }
937  else{ // Case where c is not up to date
938  assert(c_cf != 0x0) ;
939  if (vi.c_cf != 0x0) {
940  *c_cf -= *(vi.c_cf) ; // -= Mtbl_cf
941  delete c ;
942  c = 0x0 ;
943  }
944  else {
945  assert(vi.c != 0x0) ;
946  coef_i() ;
947  *c -= *(vi.c) ; // -= Mtbl
948  delete c_cf ;
949  c_cf = 0x0 ;
950  }
951  }
952 
953  // Menage (a ne faire qu'a la fin seulement)
954  del_deriv() ;
955 
956 }
957 
958 void Valeur::operator*=(const Valeur & vi) {
959 
960  // Protection
961  assert(mg == vi.get_mg()) ; // meme grille
962  assert(etat != ETATNONDEF) ; // etat defini
963  assert(vi.get_etat() != ETATNONDEF) ; // etat defini
964 
965  // Cas particuliers
966  if (etat == ETATZERO) {
967  return ;
968  }
969  if (vi.get_etat() == ETATZERO) {
970  set_etat_zero() ;
971  return ;
972  }
973 
974  // Cas general
975 
976  // Calcul dans l'espace physique
977  if (c == 0x0) {
978  coef_i() ;
979  }
980  if (vi.c == 0x0) {
981  vi.coef_i() ;
982  }
983 
984  // Calcul
985  *c *= *(vi.c) ;
986 
987  // Affectation de la base :
988  base = base * vi.base ;
989 
990  // Coefficients
991  delete c_cf ;
992  c_cf = 0x0 ;
993 
994  // Menage (a ne faire qu'a la fin seulement)
995  del_deriv() ;
996 
997 }
998 
999  //-----------------------------------//
1000  // Multiplication without aliasing //
1001  //-----------------------------------//
1002 
1003 Valeur operator%(const Valeur& t1, const Valeur& t2)
1004 {
1005  // Protections
1006  assert(t1.get_etat() != ETATNONDEF) ;
1007  assert(t2.get_etat() != ETATNONDEF) ;
1008  assert(t1.get_mg() == t2.get_mg()) ;
1009 
1010  // Cas particuliers
1011  if (t1.get_etat() == ETATZERO) {
1012  return t1 ;
1013  }
1014  if (t2.get_etat() == ETATZERO) {
1015  return t2 ;
1016  }
1017  assert(t1.get_etat() == ETATQCQ) ; // sinon...
1018  assert(t2.get_etat() == ETATQCQ) ; // sinon...
1019 
1020  // Grid
1021  const Mg3d& mg = *(t1.get_mg()) ;
1022 
1023  // Grid with twice the number of points in each dimension:
1024  const Mg3d& mg2 = *(mg.get_twice()) ;
1025 
1026  // The coefficients are required
1027  if (t1.c_cf == 0x0) {
1028  t1.coef() ;
1029  }
1030  if (t2.c_cf == 0x0) {
1031  t2.coef() ;
1032  }
1033 
1034  const Mtbl_cf& c1 = *(t1.c_cf) ;
1035  const Mtbl_cf& c2 = *(t2.c_cf) ;
1036 
1037  assert( c1.get_etat() == ETATQCQ ) ;
1038  assert( c2.get_etat() == ETATQCQ ) ;
1039 
1040  // The number of coefficients is multiplied by 2 and the additionnal
1041  // coefficients are set to zero
1042  // -----------------------------------------------------------------
1043 
1044  Mtbl_cf cc1( mg2, c1.base ) ;
1045  Mtbl_cf cc2( mg2, c2.base ) ;
1046 
1047  cc1.set_etat_qcq() ;
1048  cc2.set_etat_qcq() ;
1049 
1050  for (int l=0; l<mg.get_nzone(); l++) {
1051 
1052  int nr = mg.get_nr(l) ;
1053  int nt = mg.get_nt(l) ;
1054  int np = mg.get_np(l) ;
1055  int nr2 = mg2.get_nr(l) ;
1056  int nt2 = mg2.get_nt(l) ;
1057  int np2 = mg2.get_np(l) ;
1058 
1059  // Copy of the coefficients of t1
1060  // ------------------------------
1061 
1062  if ( c1.t[l]->get_etat() == ETATZERO ) {
1063  cc1.t[l]->set_etat_zero() ;
1064  }
1065  else {
1066 
1067  assert( c1.t[l]->get_etat() == ETATQCQ ) ;
1068  cc1.t[l]->set_etat_qcq() ;
1069 
1070  // Copy of the coefficients of t1
1071  for (int k=0; k<np+1; k++) {
1072  for (int j=0; j<nt; j++) {
1073  for (int i=0; i<nr; i++) {
1074  cc1.t[l]->set(k, j, i) = (*(c1.t[l]))(k, j, i) ;
1075  }
1076  }
1077  }
1078 
1079  // The extra phi coefficients are set to zero
1080  for (int k=np+1; k<np2+2; k++) {
1081  for (int j=0; j<nt2; j++) {
1082  for (int i=0; i<nr2; i++) {
1083  cc1.t[l]->set(k, j, i) = 0 ;
1084  }
1085  }
1086  }
1087 
1088  // The extra theta coefficients are set to zero
1089  for (int k=0; k<np+1; k++) {
1090  for (int j=nt; j<nt2; j++) {
1091  for (int i=0; i<nr2; i++) {
1092  cc1.t[l]->set(k, j, i) = 0 ;
1093  }
1094  }
1095  }
1096 
1097  // The extra r coefficients are set to zero
1098  for (int k=0; k<np+1; k++) {
1099  for (int j=0; j<nt; j++) {
1100  for (int i=nr; i<nr2; i++) {
1101  cc1.t[l]->set(k, j, i) = 0 ;
1102  }
1103  }
1104  }
1105 
1106  }
1107 
1108  // Copy of the coefficients of t2
1109  // ------------------------------
1110 
1111  if ( c2.t[l]->get_etat() == ETATZERO ) {
1112  cc2.t[l]->set_etat_zero() ;
1113  }
1114  else {
1115 
1116  assert( c2.t[l]->get_etat() == ETATQCQ ) ;
1117  cc2.t[l]->set_etat_qcq() ;
1118 
1119  // Copy of the coefficients of t2
1120  for (int k=0; k<np+1; k++) {
1121  for (int j=0; j<nt; j++) {
1122  for (int i=0; i<nr; i++) {
1123  cc2.t[l]->set(k, j, i) = (*(c2.t[l]))(k, j, i) ;
1124  }
1125  }
1126  }
1127 
1128  // The extra phi coefficients are set to zero
1129  for (int k=np+1; k<np2+2; k++) {
1130  for (int j=0; j<nt2; j++) {
1131  for (int i=0; i<nr2; i++) {
1132  cc2.t[l]->set(k, j, i) = 0 ;
1133  }
1134  }
1135  }
1136 
1137  // The extra theta coefficients are set to zero
1138  for (int k=0; k<np+1; k++) {
1139  for (int j=nt; j<nt2; j++) {
1140  for (int i=0; i<nr2; i++) {
1141  cc2.t[l]->set(k, j, i) = 0 ;
1142  }
1143  }
1144  }
1145 
1146  // The extra r coefficients are set to zero
1147  for (int k=0; k<np+1; k++) {
1148  for (int j=0; j<nt; j++) {
1149  for (int i=nr; i<nr2; i++) {
1150  cc2.t[l]->set(k, j, i) = 0 ;
1151  }
1152  }
1153  }
1154 
1155  }
1156 
1157  } // End of loop on the domains
1158 
1159 
1160  Valeur tt1( mg2 ) ;
1161  Valeur tt2( mg2 ) ;
1162 
1163  tt1 = cc1 ;
1164  tt2 = cc2 ;
1165 
1166 
1167  // Multiplication (in the configuration space) on the large grids
1168  // --------------------------------------------------------------
1169 
1170  tt1 = tt1 * tt2 ;
1171 
1172  // Coefficients of the result
1173  // --------------------------
1174 
1175  tt1.coef() ;
1176 
1177  tt1.ylm() ;
1178 
1179  const Mtbl_cf& cr2 = *(tt1.c_cf) ;
1180 
1181  Mtbl_cf cr(mg, cr2.base ) ;
1182 
1183  cr.set_etat_qcq() ;
1184 
1185  for (int l=0; l<mg.get_nzone(); l++) {
1186 
1187  if ( cr2.t[l]->get_etat() == ETATZERO ) {
1188 
1189  cr.t[l]->set_etat_zero() ;
1190 
1191  }
1192  else {
1193 
1194  assert( cr2.t[l]->get_etat() == ETATQCQ ) ;
1195 
1196  cr.t[l]->set_etat_qcq() ;
1197 
1198  int nr = mg.get_nr(l) ;
1199  int nt = mg.get_nt(l) ;
1200  int np = mg.get_np(l) ;
1201 
1202  // Copy of the coefficients of cr2
1203  for (int k=0; k<np+1; k++) {
1204  for (int j=0; j<nt; j++) {
1205  for (int i=0; i<nr; i++) {
1206  cr.t[l]->set(k, j, i) = (*(cr2.t[l]))(k, j, i) ;
1207  }
1208  }
1209  }
1210 
1211  // The last coefficient in phi is set to zero
1212  for (int j=0; j<nt; j++) {
1213  for (int i=0; i<nr; i++) {
1214  cr.t[l]->set(np+1, j, i) = 0 ;
1215  }
1216  }
1217 
1218  }
1219 
1220  } // End of loop on the domains
1221 
1222  Valeur resu( mg ) ;
1223 
1224  resu = cr ;
1225 
1226  resu.ylm_i() ;
1227 
1228  return resu ;
1229 }
1230 
1231  //---------------------------------------//
1232  // Multiplication with de-aliasing in r //
1233  //---------------------------------------//
1234 
1235 Valeur operator|(const Valeur& t1, const Valeur& t2)
1236 {
1237  // Protections
1238  assert(t1.get_etat() != ETATNONDEF) ;
1239  assert(t2.get_etat() != ETATNONDEF) ;
1240  assert(t1.get_mg() == t2.get_mg()) ;
1241 
1242  // Cas particuliers
1243  if (t1.get_etat() == ETATZERO) {
1244  return t1 ;
1245  }
1246  if (t2.get_etat() == ETATZERO) {
1247  return t2 ;
1248  }
1249  assert(t1.get_etat() == ETATQCQ) ; // sinon...
1250  assert(t2.get_etat() == ETATQCQ) ; // sinon...
1251 
1252  // Grid
1253  const Mg3d& mg = *(t1.get_mg()) ;
1254 
1255  // Grid with twice the number of points in each dimension:
1256  const Mg3d& mg2 = *(mg.plus_half()) ;
1257 
1258  // The coefficients are required
1259  if (t1.c_cf == 0x0) {
1260  t1.coef() ;
1261  }
1262  if (t2.c_cf == 0x0) {
1263  t2.coef() ;
1264  }
1265 
1266  const Mtbl_cf& c1 = *(t1.c_cf) ;
1267  const Mtbl_cf& c2 = *(t2.c_cf) ;
1268 
1269  assert( c1.get_etat() == ETATQCQ ) ;
1270  assert( c2.get_etat() == ETATQCQ ) ;
1271 
1272  // The number of coefficients is increased by 50% in r
1273  // and the additionnal coefficients are set to zero
1274  // ---------------------------------------------------
1275 
1276  Mtbl_cf cc1( mg2, c1.base ) ;
1277  Mtbl_cf cc2( mg2, c2.base ) ;
1278 
1279  cc1.set_etat_qcq() ;
1280  cc2.set_etat_qcq() ;
1281 
1282  for (int l=0; l<mg.get_nzone(); l++) {
1283 
1284  int nr = mg.get_nr(l) ;
1285  int nt = mg.get_nt(l) ;
1286  int np = mg.get_np(l) ;
1287  int nr2 = mg2.get_nr(l) ;
1288 
1289  // Copy of the coefficients of t1
1290  // ------------------------------
1291 
1292  if ( c1.t[l]->get_etat() == ETATZERO ) {
1293  cc1.t[l]->set_etat_zero() ;
1294  }
1295  else {
1296 
1297  assert( c1.t[l]->get_etat() == ETATQCQ ) ;
1298  cc1.t[l]->set_etat_qcq() ;
1299 
1300  // Copy of the coefficients of t1
1301  for (int k=0; k<np+1; k++) {
1302  for (int j=0; j<nt; j++) {
1303  for (int i=0; i<nr; i++) {
1304  cc1.t[l]->set(k, j, i) = (*(c1.t[l]))(k, j, i) ;
1305  }
1306  }
1307  }
1308 
1309  // The extra phi coefficient is set to zero
1310  for (int j=0; j<nt; j++) {
1311  for (int i=0; i<nr2; i++) {
1312  cc1.t[l]->set(np+1, j, i) = 0 ;
1313  }
1314  }
1315 
1316 
1317  // The extra r coefficients are set to zero
1318  for (int k=0; k<np+1; k++) {
1319  for (int j=0; j<nt; j++) {
1320  for (int i=nr; i<nr2; i++) {
1321  cc1.t[l]->set(k, j, i) = 0 ;
1322  }
1323  }
1324  }
1325 
1326  }
1327 
1328  // Copy of the coefficients of t2
1329  // ------------------------------
1330 
1331  if ( c2.t[l]->get_etat() == ETATZERO ) {
1332  cc2.t[l]->set_etat_zero() ;
1333  }
1334  else {
1335 
1336  assert( c2.t[l]->get_etat() == ETATQCQ ) ;
1337  cc2.t[l]->set_etat_qcq() ;
1338 
1339  // Copy of the coefficients of t2
1340  for (int k=0; k<np+1; k++) {
1341  for (int j=0; j<nt; j++) {
1342  for (int i=0; i<nr; i++) {
1343  cc2.t[l]->set(k, j, i) = (*(c2.t[l]))(k, j, i) ;
1344  }
1345  }
1346  }
1347 
1348  // The extra phi coefficient is set to zero
1349  for (int j=0; j<nt; j++) {
1350  for (int i=0; i<nr2; i++) {
1351  cc2.t[l]->set(np+1, j, i) = 0 ;
1352  }
1353  }
1354 
1355  // The extra r coefficients are set to zero
1356  for (int k=0; k<np+1; k++) {
1357  for (int j=0; j<nt; j++) {
1358  for (int i=nr; i<nr2; i++) {
1359  cc2.t[l]->set(k, j, i) = 0 ;
1360  }
1361  }
1362  }
1363 
1364  }
1365 
1366  } // End of loop on the domains
1367 
1368 
1369  Valeur tt1( mg2 ) ;
1370  Valeur tt2( mg2 ) ;
1371 
1372  tt1 = cc1 ;
1373  tt2 = cc2 ;
1374 
1375  // Multiplication (in the configuration space) on the large grids
1376  // --------------------------------------------------------------
1377 
1378  tt1 = tt1 * tt2 ;
1379 
1380  // Coefficients of the result
1381  // --------------------------
1382 
1383  tt1.coef() ;
1384 
1385  const Mtbl_cf& cr2 = *(tt1.c_cf) ;
1386 
1387  Mtbl_cf cr(mg, cr2.base ) ;
1388 
1389  cr.set_etat_qcq() ;
1390 
1391  for (int l=0; l<mg.get_nzone(); l++) {
1392 
1393  if ( cr2.t[l]->get_etat() == ETATZERO ) {
1394 
1395  cr.t[l]->set_etat_zero() ;
1396 
1397  }
1398  else {
1399 
1400  assert( cr2.t[l]->get_etat() == ETATQCQ ) ;
1401 
1402  cr.t[l]->set_etat_qcq() ;
1403 
1404  int nr = mg.get_nr(l) ;
1405  int nt = mg.get_nt(l) ;
1406  int np = mg.get_np(l) ;
1407 
1408  // Copy of the coefficients of cr2
1409  for (int k=0; k<np+1; k++) {
1410  for (int j=0; j<nt; j++) {
1411  for (int i=0; i<nr; i++) {
1412  cr.t[l]->set(k, j, i) = (*(cr2.t[l]))(k, j, i) ;
1413  }
1414  }
1415  }
1416 
1417  // The last coefficient in phi is set to zero
1418  for (int j=0; j<nt; j++) {
1419  for (int i=0; i<nr; i++) {
1420  cr.t[l]->set(np+1, j, i) = 0 ;
1421  }
1422  }
1423 
1424  }
1425 
1426  } // End of loop on the domains
1427 
1428  Valeur resu( mg ) ;
1429 
1430  resu = cr ;
1431 
1432  return resu ;
1433 }
1434 
1435 }
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
void operator+=(const Valeur &)
+= Valeur
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: valeur.C:689
Multi-domain array.
Definition: mtbl.h:118
Lorene prototypes.
Definition: app_hor.h:64
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
void annule_hard()
Sets the Valeur to zero in a hard way.
Definition: valeur.C:723
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Base_val operator*(const Base_val &, const Base_val &)
This operator is used when calling multiplication or division of Valeur .
void coef_i() const
Computes the physical value of *this.
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
int get_etat() const
Gives the logical state.
Definition: tbl.h:394
Cmp operator%(const Cmp &, const Cmp &)
Cmp * Cmp with desaliasing.
Definition: cmp_arithm.C:364
Cmp operator/(const Cmp &, const Cmp &)
Cmp / Cmp.
Definition: cmp_arithm.C:457
Scalar operator|(const Scalar &, const Scalar &)
Scalar * Scalar with desaliasing only in r.
int get_etat() const
Returns the logical state.
Definition: mtbl_cf.h:456
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
Definition: valeur.h:729
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
int get_etat() const
Gives the logical state.
Definition: mtbl.h:277
int get_etat() const
Returns the logical state.
Definition: valeur.h:726
const Mg3d * mg
Multi-grid Mgd3 on which this is defined.
Definition: valeur.h:292
Cmp operator+(const Cmp &)
Definition: cmp_arithm.C:104
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:299
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl_cf.C:300
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: tbl.C:347
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Multi-domain grid.
Definition: grilles.h:273
Base_val std_base_scal() const
Returns the standard spectral bases for a scalar.
void operator-=(const Valeur &)
-= Valeur
void operator*=(const Valeur &)
*= Valeur
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:186
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
Base_val base
Bases of the spectral expansions.
Definition: mtbl_cf.h:200
Cmp operator-(const Cmp &)
- Cmp
Definition: cmp_arithm.C:108
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition: valeur.h:295
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
void del_deriv()
Logical destructor of the derivatives.
Definition: valeur.C:638
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s which contain the spectral coefficients in each domain...
Definition: mtbl_cf.h:205