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
31char 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
144namespace Lorene {
146
147 // Protection
148 assert(vi.get_etat() != ETATNONDEF) ;
149
150 return vi ;
151}
152
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// ---------------
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// -------------
244Valeur 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// -------------
273Valeur operator+(const Mtbl& mi, const Valeur& t1) {
274 return t1 + mi ;
275}
276
277
278// Valeur + double
279// ---------------
280Valeur 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// ---------------
308Valeur operator+(double x, const Valeur& t1) // double + Valeur
309{
310 return t1 + x ;
311}
312
313// Valeur + int
314// ------------
315Valeur operator+(const Valeur& t1, int m) // Valeur + int
316{
317 return t1 + double(m) ;
318}
319
320// int + Valeur
321// -------------
322Valeur 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// ---------------
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// -------------
392Valeur 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// -------------
421Valeur operator-(const Mtbl& mi, const Valeur& t1) {
422 return - (t1 - mi) ;
423}
424
425// Valeur - double
426// ---------------
427Valeur 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// ---------------
455Valeur operator-(double x, const Valeur& t1) // double - Valeur
456{
457 return - (t1 - x) ;
458}
459
460// Valeur - int
461// ------------
462Valeur operator-(const Valeur& t1, int m) // Valeur - int
463{
464 return t1 - double(m) ;
465}
466
467// int - Valeur
468// -------------
469Valeur 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// ---------------
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// ---------------
522Valeur 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// ---------------
552Valeur operator*(double x, const Valeur& t1) // double * Valeur
553{
554 return t1 * x ;
555}
556
557// Valeur * int
558// ------------
559Valeur operator*(const Valeur& t1, int m) // Valeur * int
560{
561 return t1 * double(m) ;
562}
563
564// int * Valeur
565// ------------
566Valeur operator*(int m, const Valeur& t1) // int * Valeur
567{
568 return t1 * double(m) ;
569}
570
571
572// Valeur * Mtbl
573// --------------
574Valeur 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// --------------
601Valeur operator*(const Mtbl& c1, const Valeur& t1) {
602
603 return t1 * c1 ;
604
605}
606
607
608// Valeur * Coord
609// --------------
610Valeur 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// --------------
637Valeur operator*(const Coord& c1, const Valeur& t1) {
638
639 return t1 * c1 ;
640
641}
642
643 //**********//
644 // DIVISION //
645 //**********//
646
647// Valeur / Valeur
648// ---------------
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// ---------------
688Valeur 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// ---------------
720Valeur 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// ------------
753Valeur operator/(const Valeur& t1, int m) {
754 return t1 / double(m) ;
755}
756
757// int / Valeur
758// ------------
759Valeur operator/(int m, const Valeur& t1) {
760 return double(m) / t1 ;
761}
762
763// Valeur / Mtbl
764// ---------------
765Valeur 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// ---------------
796Valeur 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
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
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
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
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
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}
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Multi-domain grid.
Definition grilles.h:273
const Mg3d * get_twice() const
Returns the pointer on the grid which has twice the number of points in each dimension (for desaliasi...
Definition mg3d.C:539
const Mg3d * plus_half() const
Returns the pointer on the grid which has 50% more points in r dimension (for desaliasing).
Definition mg3d.C:587
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Multi-domain array.
Definition mtbl.h:118
int get_etat() const
Gives the logical state.
Definition mtbl.h:277
Values and coefficients of a (real-value) function.
Definition valeur.h:287
void operator-=(const Valeur &)
-= Valeur
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
void operator+=(const Valeur &)
+= Valeur
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition valeur.C:689
Mtbl * c
Values of the function at the points of the multi-grid
Definition valeur.h:299
const Mg3d * mg
Multi-grid Mgd3 on which this is defined.
Definition valeur.h:292
void coef_i() const
Computes the physical value of *this.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
void operator*=(const Valeur &)
*= Valeur
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
void annule_hard()
Sets the Valeur to zero in a hard way.
Definition valeur.C:723
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition valeur.h:295
void del_deriv()
Logical destructor of the derivatives.
Definition valeur.C:638
Base_val operator*(const Base_val &, const Base_val &)
This operator is used when calling multiplication or division of Valeur .
Cmp operator-(const Cmp &)
- Cmp
Definition cmp_arithm.C:108
Cmp operator/(const Cmp &, const Cmp &)
Cmp / Cmp.
Definition cmp_arithm.C:457
Cmp operator%(const Cmp &, const Cmp &)
Cmp * Cmp with desaliasing.
Definition cmp_arithm.C:364
Cmp operator+(const Cmp &)
Definition cmp_arithm.C:104
Scalar operator|(const Scalar &, const Scalar &)
Scalar * Scalar with desaliasing only in r.
Lorene prototypes.
Definition app_hor.h:64