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