valeur_arithm.C
/*
* Arithmetical operations for class Valeur
*
*/
/*
* Copyright (c) 1999-2000 Jean-Alain Marck
* Copyright (c) 1999-2005 Eric Gourgoulhon
* Copyright (c) 1999-2001 Philippe Grandclement
* Copyright (c) 2004 Jerome Novak
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
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 $" ;
/*
* $Id: valeur_arithm.C,v 1.7 2014/10/13 08:53:49 j_novak Exp $
* $Log: valeur_arithm.C,v $
* Revision 1.7 2014/10/13 08:53:49 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.6 2014/10/06 15:13:22 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.5 2008/08/27 08:52:55 jl_cornou
* Added Jacobi(0,2) polynomials case
*
* Revision 1.4 2005/11/17 15:19:23 e_gourgoulhon
* Added Valeur + Mtbl and Valeur - Mtbl.
*
* Revision 1.3 2004/07/06 13:36:30 j_novak
* Added methods for desaliased product (operator |) only in r direction.
*
* Revision 1.2 2002/10/16 14:37:15 j_novak
* Reorganization of #include instructions of standard C++, in order to
* use experimental version 3 of gcc.
*
* Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
* LORENE
*
* Revision 2.20 2001/08/31 15:23:48 eri * operator% : traitement du cas ou le Tbl est zero dans une zone.
*
* Revision 2.19 2001/05/28 12:42:27 eric
* Passage en ylm pour le desaliasing dans operator%.
*
* Revision 2.18 2001/05/26 14:50:21 eric
* Ajout de l'operator% : produit de deux Valeur avec desaliasage.
*
* Revision 2.17 2000/01/14 14:42:47 eric
* Valeur * double : operation effectue dans l'espace des coefficients si
* la Valeur n'est pas a jour ds l'espace des config.
* Valeur / double : idem.
* += Valeur : idem.
* -= Valeur : idem.
*
* Pour += Valeur et -= Valeur les schemas sont desormais calques sur
* Valeur + Valeur et Valeur - Valeur.
*
* Revision 2.16 1999/12/10 16:33:36 eric
* Dans l'arithmetique membre (+=, -=, *=), on n'appelle desormais
* del_deriv() que tout a la fin.
*
* Revision 2.15 1999/11/30 14:12:54 phil
* gestion de base dans operator/ (double,Vlaeur)
*
* Revision 2.14 1999/11/30 12:42:10 eric
* Le membre base est desormais un objet de type Base_val et non plus
* un pointeur vers une Base_val.
*
* Revision 2.13 1999/11/29 13:28:06 eric
* *** empty log message ***
*
* Revision 2.12 1999/11/29 10:20:50 eric
* Ajout de Valeur/Mtbl et Mtbl / Valeur.
*
* Revision 2.11 1999/11/29 10:06:05 eric
* Ajout de Valeur*Mtbl et Mtbl*Valeur
*
* Revision 2.10 1999/10/26 14:40:29 phil
* On gere les bases pour *, *=, /= et /
*
* Revision 2.9 1999/10/20 15:31:28 eric
* Ajout de la plupart des fonctions arithmetiques.
*
* Revision 2.8 1999/10/18 15:07:47 eric
* La fonction membre annule() est rebaptisee annule_hard().
*
* Revision 2.7 1999/10/13 15:50:56 eric
* Depoussierage.
*
* Revision 2.6 1999/09/15 10:02:26 phil
* gestion de la base dans Valeur operator (double, const Valeur &)
*
* Revision 2.5 1999/09/14 17:18:36 phil
* aout de Valeur operator* (double, const Valeur&)
*
* Revision 2.4 1999/09/13 14:52:55 phil
* *** empty log message ***
*
* Revision 2.3 1999/04/09 12:38:05 phil
* *** empty log message ***
*
* Revision 2.2 1999/04/09 12:26:59 phil
* Ajout de valeur * coord
*
* Revision 2.1 1999/02/22 15:49:28 hyc
* *** empty log message ***
*
*
* $Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_arithm.C,v 1.7 2014/10/13 08:53:49 j_novak Exp $
*
*/
// Fichiers include
// ----------------
#include <cmath>
#include <cassert>
#include <cstdlib>
#include "mtbl.h"
#include "valeur.h"
#include "coord.h"
//********************//
// OPERATEURS UNAIRES //
//********************//
namespace Lorene {
Valeur operator+(const Valeur & vi) {
// Protection
assert(vi.get_etat() != ETATNONDEF) ;
return vi ;
}
Valeur operator-(const Valeur & vi) {
// Protection
assert(vi.get_etat() != ETATNONDEF) ;
// Cas particulier
if (vi.get_etat() == ETATZERO) {
return vi ;
}
// Cas general
assert(vi.get_etat() == ETATQCQ) ; // sinon...
Valeur resu(vi.get_mg()) ; // Valeur resultat
if (vi.c != 0x0) {
resu = - *(vi.c) ;
resu.base = vi.base ; // N'oublions pas la base...
}
else{
assert(vi.c_cf != 0x0) ;
resu = - *(vi.c_cf) ; // Dans ce cas la base est prise en
// charge par l'operator=(const Mtbl_cf&).
}
// Termine
return resu ;
}
//**********//
// ADDITION //
//**********//
// Valeur + Valeur
// ---------------
Valeur operator+(const Valeur& t1, const Valeur& t2)
{
// Protections
assert(t1.get_etat() != ETATNONDEF) ;
assert(t2.get_etat() != ETATNONDEF) ;
assert(t1.get_mg() == t2.get_mg()) ;
// Cas particuliers
if (t1.get_etat() == ETATZERO) {
return t2 ;
}
if (t2.get_etat() == ETATZERO) {
return t1 ;
}
assert(t1.get_etat() == ETATQCQ) ; // sinon...
assert(t2.get_etat() == ETATQCQ) ; // sinon...
Valeur resu(t1.get_mg()) ;
// On privelegie l'addition dans l'espace des configurations:
// ----------------------------------------------------------
if (t1.c != 0x0) {
if (t2.c != 0x0) {
resu = *(t1.c) + *(t2.c) ; // Addition des Mtbl
resu.base = t1.base ;
}
else {
assert(t2.c_cf != 0x0) ;
if (t1.c_cf != 0x0) {
resu = *(t1.c_cf) + *(t2.c_cf) ; // Addition des Mtbl_cf
}
else {
t2.coef_i() ;
resu = *(t1.c) + *(t2.c) ; // Addition des Mtbl
resu.base = t1.base ;
}
}
}
else{ // Cas ou t1.c n'est pas a jour
assert(t1.c_cf != 0x0) ;
if (t2.c_cf != 0x0) {
resu = *(t1.c_cf) + *(t2.c_cf) ; // Addition des Mtbl_cf
}
else {
assert(t2.c != 0x0) ;
t1.coef_i() ;
resu = *(t1.c) + *(t2.c) ; // Addition des Mtbl
resu.base = t1.base ;
}
}
return resu ;
}
// Valeur + Mtbl
// -------------
Valeur operator+(const Valeur& t1, const Mtbl& mi)
{
// Protections
assert(t1.get_etat() != ETATNONDEF) ;
assert(mi.get_etat() != ETATNONDEF) ;
// Cas particuliers
if (mi.get_etat() == ETATZERO) {
return t1 ;
}
Valeur resu(t1) ;
if (t1.get_etat() == ETATZERO) {
resu.set_etat_c_qcq() ;
*(resu.c) = mi ;
}
else{
assert(resu.get_etat() == ETATQCQ) ; // sinon ...
resu.coef_i() ; // l'addition se fait dans l'espace des configurations
resu.set_etat_c_qcq() ;
*(resu.c) += mi ;
}
return resu ;
}
// Mtbl + Valeur
// -------------
Valeur operator+(const Mtbl& mi, const Valeur& t1) {
return t1 + mi ;
}
// Valeur + double
// ---------------
Valeur operator+(const Valeur& t1, double x)
{
// Protections
assert(t1.get_etat() != ETATNONDEF) ;
// Cas particuliers
if (x == double(0)) {
return t1 ;
}
Valeur resu(t1) ;
if (t1.get_etat() == ETATZERO) {
resu.set_etat_c_qcq() ;
*(resu.c) = x ;
}
else{
assert(resu.get_etat() == ETATQCQ) ; // sinon ...
resu.coef_i() ; // l'addition se fait dans l'espace des configurations
resu.set_etat_c_qcq() ;
*(resu.c) = *(resu.c) + x ;
}
return resu ;
}
// double + Valeur
// ---------------
Valeur operator+(double x, const Valeur& t1) // double + Valeur
{
return t1 + x ;
}
// Valeur + int
// ------------
Valeur operator+(const Valeur& t1, int m) // Valeur + int
{
return t1 + double(m) ;
}
// int + Valeur
// -------------
Valeur operator+(int m, const Valeur& t1) // int + Valeur
{
return t1 + double(m) ;
}
//**************//
// SOUSTRACTION //
//**************//
// Valeur - Valeur
// ---------------
Valeur operator-(const Valeur& t1, const Valeur& t2)
{
// Protections
assert(t1.get_etat() != ETATNONDEF) ;
assert(t2.get_etat() != ETATNONDEF) ;
assert(t1.get_mg() == t2.get_mg()) ;
// Cas particuliers
if (t1.get_etat() == ETATZERO) {
return - t2 ;
}
if (t2.get_etat() == ETATZERO) {
return t1 ;
}
assert(t1.get_etat() == ETATQCQ) ; // sinon...
assert(t2.get_etat() == ETATQCQ) ; // sinon...
Valeur resu(t1.get_mg()) ;
// On privelegie la soustraction dans l'espace des configurations:
// ---------------------------------------------------------------
if (t1.c != 0x0) {
if (t2.c != 0x0) {
resu = *(t1.c) - *(t2.c) ; // Soustraction des Mtbl
resu.base = t1.base ;
}
else {
assert(t2.c_cf != 0x0) ;
if (t1.c_cf != 0x0) {
resu = *(t1.c_cf) - *(t2.c_cf) ; // Soustraction des Mtbl_cf
}
else {
t2.coef_i() ;
resu = *(t1.c) - *(t2.c) ; // Soustraction des Mtbl
resu.base = t1.base ;
}
}
}
else{ // Cas ou t1.c n'est pas a jour
assert(t1.c_cf != 0x0) ;
if (t2.c_cf != 0x0) {
resu = *(t1.c_cf) - *(t2.c_cf) ; // Soustraction des Mtbl_cf
}
else {
assert(t2.c != 0x0) ;
t1.coef_i() ;
resu = *(t1.c) - *(t2.c) ; // Soustraction des Mtbl
resu.base = t1.base ;
}
}
return resu ;
}
// Valeur - Mtbl
// -------------
Valeur operator-(const Valeur& t1, const Mtbl& mi)
{
// Protections
assert(t1.get_etat() != ETATNONDEF) ;
assert(mi.get_etat() != ETATNONDEF) ;
// Cas particuliers
if (mi.get_etat() == ETATZERO) {
return t1 ;
}
Valeur resu(t1) ;
if (t1.get_etat() == ETATZERO) {
resu.set_etat_c_qcq() ;
*(resu.c) = - mi ;
}
else{
assert(resu.get_etat() == ETATQCQ) ; // sinon ...
resu.coef_i() ; // substraction in configuration space
resu.set_etat_c_qcq() ;
*(resu.c) -= mi ;
}
return resu ;
}
// Mtbl - Valeur
// -------------
Valeur operator-(const Mtbl& mi, const Valeur& t1) {
return - (t1 - mi) ;
}
// Valeur - double
// ---------------
Valeur operator-(const Valeur& t1, double x)
{
// Protections
assert(t1.get_etat() != ETATNONDEF) ;
// Cas particuliers
if (x == double(0)) {
return t1 ;
}
Valeur resu(t1) ;
if (t1.get_etat() == ETATZERO) {
resu.set_etat_c_qcq() ;
*(resu.c) = - x ;
}
else{
assert(resu.get_etat() == ETATQCQ) ; // sinon ...
resu.coef_i() ; // l'addition se fait dans l'espace des configurations
resu.set_etat_c_qcq() ;
*(resu.c) = *(resu.c) - x ;
}
return resu ;
}
// double - Valeur
// ---------------
Valeur operator-(double x, const Valeur& t1) // double - Valeur
{
return - (t1 - x) ;
}
// Valeur - int
// ------------
Valeur operator-(const Valeur& t1, int m) // Valeur - int
{
return t1 - double(m) ;
}
// int - Valeur
// -------------
Valeur operator-(int m, const Valeur& t1) // int - Valeur
{
return double(m) - t1 ;
}
//****************//
// MULTIPLICATION //
//****************//
// Valeur * Valeur
// ---------------
Valeur operator*(const Valeur& t1, const Valeur& t2)
{
// Protections
assert(t1.get_etat() != ETATNONDEF) ;
assert(t2.get_etat() != ETATNONDEF) ;
assert(t1.get_mg() == t2.get_mg()) ;
// Cas particuliers
if (t1.get_etat() == ETATZERO) {
return t1 ;
}
if (t2.get_etat() == ETATZERO) {
return t2 ;
}
assert(t1.get_etat() == ETATQCQ) ; // sinon...
assert(t2.get_etat() == ETATQCQ) ; // sinon...
Valeur resu(t1.get_mg()) ;
// La multiplication est faite dans l'espace des configurations:
if (t1.c == 0x0) {
t1.coef_i() ;
}
if (t2.c == 0x0) {
t2.coef_i() ;
}
resu = (*(t1.c)) * (*(t2.c)) ; // Multiplication des Mtbl
// affectation de la base :
resu.base = t1.base * t2.base;
return resu ;
}
// Valeur * double
// ---------------
Valeur operator*(const Valeur& c1, double a) {
// Protection
assert(c1.get_etat() != ETATNONDEF) ;
// Cas particulier
if ((c1.get_etat() == ETATZERO) || ( a == double(1) )) {
return c1 ;
}
// Cas general
assert(c1.get_etat() == ETATQCQ) ; // sinon...
Valeur result(c1.get_mg()) ;
if (c1.c != 0x0) {
result = *(c1.c) * a ; // Mtbl * double
result.base = c1.base ; // in this case, result.base must be set
// by hand
}
else {
assert(c1.c_cf != 0x0) ;
result = *(c1.c_cf) * a ; // Mtbl_cf * double
}
return result ;
}
// double * Valeur
// ---------------
Valeur operator*(double x, const Valeur& t1) // double * Valeur
{
return t1 * x ;
}
// Valeur * int
// ------------
Valeur operator*(const Valeur& t1, int m) // Valeur * int
{
return t1 * double(m) ;
}
// int * Valeur
// ------------
Valeur operator*(int m, const Valeur& t1) // int * Valeur
{
return t1 * double(m) ;
}
// Valeur * Mtbl
// --------------
Valeur operator*(const Valeur & c1, const Mtbl& c2) {
// Protection
assert(c1.get_etat() != ETATNONDEF) ;
// Cas particulier
if (c1.get_etat() == ETATZERO) {
return c1 ;
}
// Cas general
assert(c1.get_etat() == ETATQCQ) ;
Valeur result(c1.get_mg()) ;
// La multiplication est faite dans l'espace des configurations:
if (c1.c == 0x0) {
c1.coef_i() ;
}
result = *(c1.c) * c2 ; // Mtbl * Mtbl
return result ;
}
// Mtbl * Valeur
// --------------
Valeur operator*(const Mtbl& c1, const Valeur& t1) {
return t1 * c1 ;
}
// Valeur * Coord
// --------------
Valeur operator*(const Valeur & c1, const Coord& c2) {
// Protection
assert(c1.get_etat() != ETATNONDEF) ;
// Cas particulier
if (c1.get_etat() == ETATZERO) {
return c1 ;
}
// Cas general
assert(c1.get_etat() == ETATQCQ) ;
Valeur result(c1.get_mg()) ;
// La multiplication est faite dans l'espace des configurations:
if (c1.c == 0x0) {
c1.coef_i() ;
}
result = *(c1.c) * c2 ; // Mtbl * Coord
return result ;
}
// Coord * Valeur
// --------------
Valeur operator*(const Coord& c1, const Valeur& t1) {
return t1 * c1 ;
}
//**********//
// DIVISION //
//**********//
// Valeur / Valeur
// ---------------
Valeur operator/(const Valeur& t1, const Valeur& t2)
{
// Protections
assert(t1.get_etat() != ETATNONDEF) ;
assert(t2.get_etat() != ETATNONDEF) ;
assert(t1.get_mg() == t2.get_mg()) ;
// Cas particuliers
if (t2.get_etat() == ETATZERO) {
cout << "Division by 0 in Valeur / Valeur !" << endl ;
abort() ;
}
if (t1.get_etat() == ETATZERO) {
return t1 ;
}
assert(t1.get_etat() == ETATQCQ) ; // sinon...
assert(t2.get_etat() == ETATQCQ) ; // sinon...
Valeur resu(t1.get_mg()) ;
// La division est faite dans l'espace des configurations:
if (t1.c == 0x0) {
t1.coef_i() ;
}
if (t2.c == 0x0) {
t2.coef_i() ;
}
resu = (*(t1.c)) / (*(t2.c)) ; // Division des Mtbl
// affectation de la base :
resu.base = t1.base * t2.base;
return resu ;
}
// Valeur / double
// ---------------
Valeur operator/(const Valeur& t1, double x)
{
// Protections
assert(t1.get_etat() != ETATNONDEF) ;
// Cas particuliers
if ( x == double(0) ) {
cout << "Division by 0 in Valeur / double !" << endl ;
abort() ;
}
if ((t1.get_etat() == ETATZERO) || ( x == double(1) )) {
return t1 ;
}
assert(t1.get_etat() == ETATQCQ) ; // sinon...
Valeur resu(t1.get_mg()) ;
if (t1.c != 0x0) {
resu = *(t1.c) / x ; // Mtbl / double
resu.base = t1.base ; // in this case, resu.base must be set by hand
}
else {
assert(t1.c_cf != 0x0) ;
resu = *(t1.c_cf) / x ; // Mtbl_cf * double
}
return resu ;
}
// double / Valeur
// ---------------
Valeur operator/(double x, const Valeur & c2) {
// Protection
assert(c2.get_etat() != ETATNONDEF) ;
// Cas particuliers
if (c2.get_etat() == ETATZERO) {
cout << "Division by 0 in double / Valeur !" << endl ;
abort() ;
}
// Cas general
assert(c2.get_etat() == ETATQCQ) ; // sinon...
// Il faut les valeurs physiques de c2
if (c2.c == 0x0) {
c2.coef_i() ;
}
// Le resultat
Valeur r(c2.get_mg()) ;
r.set_etat_c_qcq() ;
*(r.c) = x / *(c2.c) ;
// affectation de la base :
r.base = c2.get_mg()->std_base_scal() * c2.base ;
// Termine
return r ;
}
// Valeur / int
// ------------
Valeur operator/(const Valeur& t1, int m) {
return t1 / double(m) ;
}
// int / Valeur
// ------------
Valeur operator/(int m, const Valeur& t1) {
return double(m) / t1 ;
}
// Valeur / Mtbl
// ---------------
Valeur operator/(const Valeur& t1, const Mtbl& m2)
{
// Protections
assert(t1.get_etat() != ETATNONDEF) ;
assert(m2.get_etat() != ETATNONDEF) ;
// Cas particuliers
if ( m2.get_etat() == ETATZERO ) {
cout << "Division by 0 in Valeur / Mtbl !" << endl ;
abort() ;
}
if (t1.get_etat() == ETATZERO) {
return t1 ;
}
assert(t1.get_etat() == ETATQCQ) ; // sinon...
Valeur resu(t1.get_mg()) ;
// La division est faite dans l'espace des configurations:
if (t1.c == 0x0) {
t1.coef_i() ;
}
resu = (*(t1.c)) / m2 ; // Division Mtbl / Mtbl
return resu ;
}
// Mtbl / Valeur
// ---------------
Valeur operator/(const Mtbl& m1, const Valeur & c2) {
// Protection
assert(m1.get_etat() != ETATNONDEF) ;
assert(c2.get_etat() != ETATNONDEF) ;
// Cas particuliers
if (c2.get_etat() == ETATZERO) {
cout << "Division by 0 in Mtbl / Valeur !" << endl ;
abort() ;
}
// Cas general
assert(c2.get_etat() == ETATQCQ) ; // sinon...
// Le resultat
Valeur resu(c2.get_mg()) ;
// Il faut les valeurs physiques de c2
if (c2.c == 0x0) {
c2.coef_i() ;
}
resu = m1 / (*(c2.c)) ; // Division Mtbl / Mtbl
// Termine
return resu ;
}
//*******************//
// operateurs +=,... //
//*******************//
void Valeur::operator+=(const Valeur & vi) {
// Protection
assert(mg == vi.get_mg()) ; // meme grille
assert(etat != ETATNONDEF) ; // etat defini
assert(vi.get_etat() != ETATNONDEF) ; // etat defini
// Cas particulier
if (vi.get_etat() == ETATZERO) {
return ;
}
// Cas general
// Cas de l'etat ZERO
if (etat == ETATZERO) {
annule_hard() ;
}
if (c != 0x0) {
if (vi.c != 0x0) {
*c += *(vi.c) ; // += Mtbl
delete c_cf ;
c_cf = 0x0 ;
}
else {
assert(vi.c_cf != 0x0) ;
if (c_cf != 0x0) {
*c_cf += *(vi.c_cf) ; // += Mtbl_cf
delete c ;
c = 0x0 ;
}
else {
vi.coef_i() ;
*c += *(vi.c) ; // += Mtbl
delete c_cf ;
c_cf = 0x0 ;
}
}
}
else{ // Case where c is not up to date
assert(c_cf != 0x0) ;
if (vi.c_cf != 0x0) {
*c_cf += *(vi.c_cf) ; // += Mtbl_cf
delete c ;
c = 0x0 ;
}
else {
assert(vi.c != 0x0) ;
coef_i() ;
*c += *(vi.c) ; // += Mtbl
delete c_cf ;
c_cf = 0x0 ;
}
}
// Menage (a ne faire qu'a la fin seulement)
del_deriv() ;
}
void Valeur::operator-=(const Valeur & vi) {
// Protection
assert(mg == vi.get_mg()) ; // meme grille
assert(etat != ETATNONDEF) ; // etat defini
assert(vi.get_etat() != ETATNONDEF) ; // etat defini
// Cas particulier
if (vi.get_etat() == ETATZERO) {
return ;
}
// Cas general
// Cas de l'etat ZERO
if (etat == ETATZERO) {
annule_hard() ;
}
if (c != 0x0) {
if (vi.c != 0x0) {
*c -= *(vi.c) ; // -= Mtbl
delete c_cf ;
c_cf = 0x0 ;
}
else {
assert(vi.c_cf != 0x0) ;
if (c_cf != 0x0) {
*c_cf -= *(vi.c_cf) ; // -= Mtbl_cf
delete c ;
c = 0x0 ;
}
else {
vi.coef_i() ;
*c -= *(vi.c) ; // -= Mtbl
delete c_cf ;
c_cf = 0x0 ;
}
}
}
else{ // Case where c is not up to date
assert(c_cf != 0x0) ;
if (vi.c_cf != 0x0) {
*c_cf -= *(vi.c_cf) ; // -= Mtbl_cf
delete c ;
c = 0x0 ;
}
else {
assert(vi.c != 0x0) ;
coef_i() ;
*c -= *(vi.c) ; // -= Mtbl
delete c_cf ;
c_cf = 0x0 ;
}
}
// Menage (a ne faire qu'a la fin seulement)
del_deriv() ;
}
void Valeur::operator*=(const Valeur & vi) {
// Protection
assert(mg == vi.get_mg()) ; // meme grille
assert(etat != ETATNONDEF) ; // etat defini
assert(vi.get_etat() != ETATNONDEF) ; // etat defini
// Cas particuliers
if (etat == ETATZERO) {
return ;
}
if (vi.get_etat() == ETATZERO) {
set_etat_zero() ;
return ;
}
// Cas general
// Calcul dans l'espace physique
if (c == 0x0) {
coef_i() ;
}
if (vi.c == 0x0) {
vi.coef_i() ;
}
// Calcul
*c *= *(vi.c) ;
// Affectation de la base :
base = base * vi.base ;
// Coefficients
delete c_cf ;
c_cf = 0x0 ;
// Menage (a ne faire qu'a la fin seulement)
del_deriv() ;
}
//-----------------------------------//
// Multiplication without aliasing //
//-----------------------------------//
Valeur operator%(const Valeur& t1, const Valeur& t2)
{
// Protections
assert(t1.get_etat() != ETATNONDEF) ;
assert(t2.get_etat() != ETATNONDEF) ;
assert(t1.get_mg() == t2.get_mg()) ;
// Cas particuliers
if (t1.get_etat() == ETATZERO) {
return t1 ;
}
if (t2.get_etat() == ETATZERO) {
return t2 ;
}
assert(t1.get_etat() == ETATQCQ) ; // sinon...
assert(t2.get_etat() == ETATQCQ) ; // sinon...
// Grid
const Mg3d& mg = *(t1.get_mg()) ;
// Grid with twice the number of points in each dimension:
const Mg3d& mg2 = *(mg.get_twice()) ;
// The coefficients are required
if (t1.c_cf == 0x0) {
t1.coef() ;
}
if (t2.c_cf == 0x0) {
t2.coef() ;
}
const Mtbl_cf& c1 = *(t1.c_cf) ;
const Mtbl_cf& c2 = *(t2.c_cf) ;
assert( c1.get_etat() == ETATQCQ ) ;
assert( c2.get_etat() == ETATQCQ ) ;
// The number of coefficients is multiplied by 2 and the additionnal
// coefficients are set to zero
// -----------------------------------------------------------------
Mtbl_cf cc1( mg2, c1.base ) ;
Mtbl_cf cc2( mg2, c2.base ) ;
cc1.set_etat_qcq() ;
cc2.set_etat_qcq() ;
for (int l=0; l<mg.get_nzone(); l++) {
int nr = mg.get_nr(l) ;
int nt = mg.get_nt(l) ;
int np = mg.get_np(l) ;
int nr2 = mg2.get_nr(l) ;
int nt2 = mg2.get_nt(l) ;
int np2 = mg2.get_np(l) ;
// Copy of the coefficients of t1
// ------------------------------
if ( c1.t[l]->get_etat() == ETATZERO ) {
cc1.t[l]->set_etat_zero() ;
}
else {
assert( c1.t[l]->get_etat() == ETATQCQ ) ;
cc1.t[l]->set_etat_qcq() ;
// Copy of the coefficients of t1
for (int k=0; k<np+1; k++) {
for (int j=0; j<nt; j++) {
for (int i=0; i<nr; i++) {
cc1.t[l]->set(k, j, i) = (*(c1.t[l]))(k, j, i) ;
}
}
}
// The extra phi coefficients are set to zero
for (int k=np+1; k<np2+2; k++) {
for (int j=0; j<nt2; j++) {
for (int i=0; i<nr2; i++) {
cc1.t[l]->set(k, j, i) = 0 ;
}
}
}
// The extra theta coefficients are set to zero
for (int k=0; k<np+1; k++) {
for (int j=nt; j<nt2; j++) {
for (int i=0; i<nr2; i++) {
cc1.t[l]->set(k, j, i) = 0 ;
}
}
}
// The extra r coefficients are set to zero
for (int k=0; k<np+1; k++) {
for (int j=0; j<nt; j++) {
for (int i=nr; i<nr2; i++) {
cc1.t[l]->set(k, j, i) = 0 ;
}
}
}
}
// Copy of the coefficients of t2
// ------------------------------
if ( c2.t[l]->get_etat() == ETATZERO ) {
cc2.t[l]->set_etat_zero() ;
}
else {
assert( c2.t[l]->get_etat() == ETATQCQ ) ;
cc2.t[l]->set_etat_qcq() ;
// Copy of the coefficients of t2
for (int k=0; k<np+1; k++) {
for (int j=0; j<nt; j++) {
for (int i=0; i<nr; i++) {
cc2.t[l]->set(k, j, i) = (*(c2.t[l]))(k, j, i) ;
}
}
}
// The extra phi coefficients are set to zero
for (int k=np+1; k<np2+2; k++) {
for (int j=0; j<nt2; j++) {
for (int i=0; i<nr2; i++) {
cc2.t[l]->set(k, j, i) = 0 ;
}
}
}
// The extra theta coefficients are set to zero
for (int k=0; k<np+1; k++) {
for (int j=nt; j<nt2; j++) {
for (int i=0; i<nr2; i++) {
cc2.t[l]->set(k, j, i) = 0 ;
}
}
}
// The extra r coefficients are set to zero
for (int k=0; k<np+1; k++) {
for (int j=0; j<nt; j++) {
for (int i=nr; i<nr2; i++) {
cc2.t[l]->set(k, j, i) = 0 ;
}
}
}
}
} // End of loop on the domains
Valeur tt1( mg2 ) ;
Valeur tt2( mg2 ) ;
tt1 = cc1 ;
tt2 = cc2 ;
// Multiplication (in the configuration space) on the large grids
// --------------------------------------------------------------
tt1 = tt1 * tt2 ;
// Coefficients of the result
// --------------------------
tt1.coef() ;
tt1.ylm() ;
const Mtbl_cf& cr2 = *(tt1.c_cf) ;
Mtbl_cf cr(mg, cr2.base ) ;
cr.set_etat_qcq() ;
for (int l=0; l<mg.get_nzone(); l++) {
if ( cr2.t[l]->get_etat() == ETATZERO ) {
cr.t[l]->set_etat_zero() ;
}
else {
assert( cr2.t[l]->get_etat() == ETATQCQ ) ;
cr.t[l]->set_etat_qcq() ;
int nr = mg.get_nr(l) ;
int nt = mg.get_nt(l) ;
int np = mg.get_np(l) ;
// Copy of the coefficients of cr2
for (int k=0; k<np+1; k++) {
for (int j=0; j<nt; j++) {
for (int i=0; i<nr; i++) {
cr.t[l]->set(k, j, i) = (*(cr2.t[l]))(k, j, i) ;
}
}
}
// The last coefficient in phi is set to zero
for (int j=0; j<nt; j++) {
for (int i=0; i<nr; i++) {
cr.t[l]->set(np+1, j, i) = 0 ;
}
}
}
} // End of loop on the domains
Valeur resu( mg ) ;
resu = cr ;
resu.ylm_i() ;
return resu ;
}
//---------------------------------------//
// Multiplication with de-aliasing in r //
//---------------------------------------//
Valeur operator|(const Valeur& t1, const Valeur& t2)
{
// Protections
assert(t1.get_etat() != ETATNONDEF) ;
assert(t2.get_etat() != ETATNONDEF) ;
assert(t1.get_mg() == t2.get_mg()) ;
// Cas particuliers
if (t1.get_etat() == ETATZERO) {
return t1 ;
}
if (t2.get_etat() == ETATZERO) {
return t2 ;
}
assert(t1.get_etat() == ETATQCQ) ; // sinon...
assert(t2.get_etat() == ETATQCQ) ; // sinon...
// Grid
const Mg3d& mg = *(t1.get_mg()) ;
// Grid with twice the number of points in each dimension:
const Mg3d& mg2 = *(mg.plus_half()) ;
// The coefficients are required
if (t1.c_cf == 0x0) {
t1.coef() ;
}
if (t2.c_cf == 0x0) {
t2.coef() ;
}
const Mtbl_cf& c1 = *(t1.c_cf) ;
const Mtbl_cf& c2 = *(t2.c_cf) ;
assert( c1.get_etat() == ETATQCQ ) ;
assert( c2.get_etat() == ETATQCQ ) ;
// The number of coefficients is increased by 50% in r
// and the additionnal coefficients are set to zero
// ---------------------------------------------------
Mtbl_cf cc1( mg2, c1.base ) ;
Mtbl_cf cc2( mg2, c2.base ) ;
cc1.set_etat_qcq() ;
cc2.set_etat_qcq() ;
for (int l=0; l<mg.get_nzone(); l++) {
int nr = mg.get_nr(l) ;
int nt = mg.get_nt(l) ;
int np = mg.get_np(l) ;
int nr2 = mg2.get_nr(l) ;
// Copy of the coefficients of t1
// ------------------------------
if ( c1.t[l]->get_etat() == ETATZERO ) {
cc1.t[l]->set_etat_zero() ;
}
else {
assert( c1.t[l]->get_etat() == ETATQCQ ) ;
cc1.t[l]->set_etat_qcq() ;
// Copy of the coefficients of t1
for (int k=0; k<np+1; k++) {
for (int j=0; j<nt; j++) {
for (int i=0; i<nr; i++) {
cc1.t[l]->set(k, j, i) = (*(c1.t[l]))(k, j, i) ;
}
}
}
// The extra phi coefficient is set to zero
for (int j=0; j<nt; j++) {
for (int i=0; i<nr2; i++) {
cc1.t[l]->set(np+1, j, i) = 0 ;
}
}
// The extra r coefficients are set to zero
for (int k=0; k<np+1; k++) {
for (int j=0; j<nt; j++) {
for (int i=nr; i<nr2; i++) {
cc1.t[l]->set(k, j, i) = 0 ;
}
}
}
}
// Copy of the coefficients of t2
// ------------------------------
if ( c2.t[l]->get_etat() == ETATZERO ) {
cc2.t[l]->set_etat_zero() ;
}
else {
assert( c2.t[l]->get_etat() == ETATQCQ ) ;
cc2.t[l]->set_etat_qcq() ;
// Copy of the coefficients of t2
for (int k=0; k<np+1; k++) {
for (int j=0; j<nt; j++) {
for (int i=0; i<nr; i++) {
cc2.t[l]->set(k, j, i) = (*(c2.t[l]))(k, j, i) ;
}
}
}
// The extra phi coefficient is set to zero
for (int j=0; j<nt; j++) {
for (int i=0; i<nr2; i++) {
cc2.t[l]->set(np+1, j, i) = 0 ;
}
}
// The extra r coefficients are set to zero
for (int k=0; k<np+1; k++) {
for (int j=0; j<nt; j++) {
for (int i=nr; i<nr2; i++) {
cc2.t[l]->set(k, j, i) = 0 ;
}
}
}
}
} // End of loop on the domains
Valeur tt1( mg2 ) ;
Valeur tt2( mg2 ) ;
tt1 = cc1 ;
tt2 = cc2 ;
// Multiplication (in the configuration space) on the large grids
// --------------------------------------------------------------
tt1 = tt1 * tt2 ;
// Coefficients of the result
// --------------------------
tt1.coef() ;
const Mtbl_cf& cr2 = *(tt1.c_cf) ;
Mtbl_cf cr(mg, cr2.base ) ;
cr.set_etat_qcq() ;
for (int l=0; l<mg.get_nzone(); l++) {
if ( cr2.t[l]->get_etat() == ETATZERO ) {
cr.t[l]->set_etat_zero() ;
}
else {
assert( cr2.t[l]->get_etat() == ETATQCQ ) ;
cr.t[l]->set_etat_qcq() ;
int nr = mg.get_nr(l) ;
int nt = mg.get_nt(l) ;
int np = mg.get_np(l) ;
// Copy of the coefficients of cr2
for (int k=0; k<np+1; k++) {
for (int j=0; j<nt; j++) {
for (int i=0; i<nr; i++) {
cr.t[l]->set(k, j, i) = (*(cr2.t[l]))(k, j, i) ;
}
}
}
// The last coefficient in phi is set to zero
for (int j=0; j<nt; j++) {
for (int i=0; i<nr; i++) {
cr.t[l]->set(np+1, j, i) = 0 ;
}
}
}
} // End of loop on the domains
Valeur resu( mg ) ;
resu = cr ;
return resu ;
}
}
Lorene::Coord
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
Lorene::Mg3d::get_twice
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
Lorene::Mg3d::plus_half
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
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Lorene::Mg3d::std_base_scal
Base_val std_base_scal() const
Returns the standard spectral bases for a scalar.
Definition mg3d_std_base.C:101
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Lorene::Mtbl_cf::base
Base_val base
Bases of the spectral expansions.
Definition mtbl_cf.h:200
Lorene::Mtbl_cf::get_etat
int get_etat() const
Returns the logical state.
Definition mtbl_cf.h:456
Lorene::Mtbl_cf::t
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
Definition mtbl_cf.h:205
Lorene::Mtbl_cf::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition mtbl_cf.C:300
Lorene::Mtbl
Multi-domain array.
Definition mtbl.h:118
Lorene::Mtbl::get_etat
int get_etat() const
Gives the logical state.
Definition mtbl.h:277
Lorene::Tbl::get_etat
int get_etat() const
Gives the logical state.
Definition tbl.h:394
Lorene::Tbl::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition tbl.C:347
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Lorene::Valeur::operator-=
void operator-=(const Valeur &)
-= Valeur
Definition valeur_arithm.C:897
Lorene::Valeur::set_etat_c_qcq
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c ).
Definition valeur.C:701
Lorene::Valeur::get_etat
int get_etat() const
Returns the logical state.
Definition valeur.h:726
Lorene::Valeur::operator+=
void operator+=(const Valeur &)
+= Valeur
Definition valeur_arithm.C:834
Lorene::Valeur::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition valeur.C:689
Lorene::Valeur::Valeur
Valeur(const Mg3d &mgrid)
Constructor.
Definition valeur.C:200
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Lorene::Valeur::c
Mtbl * c
Values of the function at the points of the multi-grid.
Definition valeur.h:299
Lorene::Valeur::mg
const Mg3d * mg
Multi-grid Mgd3 on which this is defined.
Definition valeur.h:292
Lorene::Valeur::coef_i
void coef_i() const
Computes the physical value of *this.
Definition valeur_coef_i.C:140
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition valeur_coef.C:148
Lorene::Valeur::operator*=
void operator*=(const Valeur &)
*= Valeur
Definition valeur_arithm.C:958
Lorene::Valeur::get_mg
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
Definition valeur.h:729
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm().
Definition valeur_ylm_i.C:131
Lorene::Valeur::base
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
Lorene::Valeur::annule_hard
void annule_hard()
Sets the Valeur to zero in a hard way.
Definition valeur.C:723
Lorene::Valeur::etat
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition valeur.h:295
Lorene::Valeur::del_deriv
void del_deriv()
Logical destructor of the derivatives.
Definition valeur.C:638
Lorene::operator*
Base_val operator*(const Base_val &, const Base_val &)
This operator is used when calling multiplication or division of Valeur .
Definition base_val_mult.C:102
Lorene::operator-
Cmp operator-(const Cmp &)
- Cmp
Definition cmp_arithm.C:108
Lorene::operator/
Cmp operator/(const Cmp &, const Cmp &)
Cmp / Cmp.
Definition cmp_arithm.C:457
Lorene::operator%
Cmp operator%(const Cmp &, const Cmp &)
Cmp * Cmp with desaliasing.
Definition cmp_arithm.C:364
Lorene::operator+
Cmp operator+(const Cmp &)
Definition cmp_arithm.C:104
Lorene::operator|
Scalar operator|(const Scalar &, const Scalar &)
Scalar * Scalar with desaliasing only in r.
Definition scalar_arithm.C:491
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::x
Coord x
x coordinate centered on the grid
Definition map.h:726
Lorene::r
Coord r
r coordinate centered on the grid
Definition map.h:718