valeur_val_propre_1d.C
/*
* Copyright (c) 2004 Philippe Grandclement
*
* 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_val_propre_1d_C[] = "$Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_val_propre_1d.C,v 1.3 2014/10/13 08:53:51 j_novak Exp $" ;
/*
* $Id: valeur_val_propre_1d.C,v 1.3 2014/10/13 08:53:51 j_novak Exp $
* $Log: valeur_val_propre_1d.C,v $
* Revision 1.3 2014/10/13 08:53:51 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.2 2014/10/06 15:13:24 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.1 2004/08/24 09:14:52 p_grandclement
* Addition of some new operators, like Poisson in 2d... It now requieres the
* GSL library to work.
*
* Also, the way a variable change is stored by a Param_elliptic is changed and
* no longer uses Change_var but rather 2 Scalars. The codes using that feature
* will requiere some modification. (It should concern only the ones about monopoles)
*
*
* $Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_val_propre_1d.C,v 1.3 2014/10/13 08:53:51 j_novak Exp $
*
*/
// headers du C
#include <cassert>
#include <cstdlib>
// headers Lorene
#include "type_parite.h"
#include "valeur.h"
#include "matrice.h"
#include "proto.h"
namespace Lorene {
//****************************************************************
// CAS PAIR
//****************************************************************
void rotate_propre_pair (Valeur& so, bool sens) {
so.coef() ;
so.set_etat_cf_qcq() ;
static int nt_courant = 0 ;
static Matrice* passage = 0x0 ;
static Matrice* inv = 0x0 ;
int nt = so.mg->get_nt(0) ;
if (nt != nt_courant) {
// On doit calculer les matrices... Pas de la tarte...
if (passage != 0x0)
delete passage ;
if (inv != 0x0)
delete inv ;
nt_courant = nt ;
Matrice ope (nt, nt) ;
ope.set_etat_qcq() ;
for (int i=0 ; i<nt ; i++)
for (int j=0 ; j<nt ; j++)
ope.set(i,j) = 0 ;
double c_courant = 0 ;
for (int j=0 ; j<nt ; j++) {
ope.set(0, j) = 2*j ;
for (int i=1 ; i<j ; i++)
ope.set(i,j) = 4*j ;
ope.set(j,j) = c_courant ;
c_courant -= 8*j + 2 ;
}
passage = new Matrice (ope.vect_propre()) ;
passage->set_band(nt-1,0) ;
passage->set_lu() ;
// Un peu nul pour calculer l'inverse mais bon...
inv = new Matrice (nt, nt) ;
inv->set_etat_qcq() ;
Tbl auxi (nt) ;
auxi.set_etat_qcq() ;
for (int i=0 ; i<nt ; i++) {
for (int j=0 ; j<nt ; j++)
auxi.set(j) = 0 ;
auxi.set(i) = 1 ;
Tbl sortie (passage->inverse(auxi)) ;
for (int j=0 ; j<nt ; j++)
inv->set(j,i) = sortie(j) ;
}
}
// Fin du calcul des matrices...
for (int l=0 ; l<so.mg->get_nzone() ; l++)
if (so.c_cf->t[l]->get_etat() != ETATZERO)
for (int k=0 ; k<so.mg->get_np(l) ; k++)
for (int i=0 ; i<so.mg->get_nr(l) ; i++) {
Tbl coefs(nt) ;
coefs.set_etat_qcq() ;
for (int j=0 ; j<nt ; j++)
coefs.set(j) = (*so.c_cf)(l,k,j,i) ;
Tbl prod (nt) ;
prod.set_etat_qcq() ;
for (int j=0 ; j<nt ; j++)
prod.set(j) = 0 ;
if (sens) {
//calcul direct
for (int j=0 ; j<nt ; j++)
for (int jb=0 ; jb<nt ; jb++)
prod.set(j) += (*inv)(j, jb) * coefs(jb) ;
}
else {
//calcul inverse
for (int j=0 ; j<nt ; j++)
for (int jb=0 ; jb<nt ; jb++)
prod.set(j) += (*passage)(j, jb) * coefs(jb) ;
}
for (int j=0 ; j<nt ; j++)
so.c_cf->set(l,k,j,i) = prod(j) ;
}
// On met les bonnes bases :
int base_tet = so.base.get_base_t(0) ;
Base_val new_base (so.base) ;
if (sens) {
switch (base_tet) {
case T_COS_P:
new_base.set_base_t(T_CL_COS_P) ;
break ;
case T_SIN_P:
new_base.set_base_t(T_CL_SIN_P) ;
break ;
default:
cout << "Problem in rotate_propre_pair" << endl ;
abort() ;
break ;
}
}
else {
switch (base_tet) {
case T_CL_COS_P:
new_base.set_base_t(T_COS_P) ;
break ;
case T_CL_SIN_P:
new_base.set_base_t(T_SIN_P) ;
break ;
default:
cout << "Problem in rotate_propre_pair" << endl ;
abort() ;
break ;
}
}
so.c_cf->base = new_base ;
so.base = new_base ;
}
//****************************************************************
// CAS IMPAIR
//****************************************************************
void rotate_propre_impair (Valeur& so, bool sens) {
so.coef() ;
so.set_etat_cf_qcq() ;
static int nt_courant = 0 ;
static Matrice* passage = 0x0 ;
static Matrice* inv = 0x0 ;
int nt = so.mg->get_nt(0) ;
if (nt != nt_courant) {
// On doit calculer les matrices... Pas de la tarte...
if (passage != 0x0)
delete passage ;
if (inv != 0x0)
delete inv ;
nt_courant = nt ;
Matrice ope (nt, nt) ;
ope.set_etat_qcq() ;
for (int i=0 ; i<nt ; i++)
for (int j=0 ; j<nt ; j++)
ope.set(i,j) = 0 ;
double c_courant = 0 ;
for (int j=0 ; j<nt ; j++) {
for (int i=0 ; i<j ; i++)
ope.set(i,j) = 2+4*j ;
ope.set(j,j) = c_courant ;
c_courant -= 8*j + 6 ;
}
passage = new Matrice (ope.vect_propre()) ;
passage->set_band(nt-1,0) ;
passage->set_lu() ;
// Un peu nul pour calculer l'inverse mais bon...
inv = new Matrice (nt, nt) ;
inv->set_etat_qcq() ;
Tbl auxi (nt) ;
auxi.set_etat_qcq() ;
for (int i=0 ; i<nt ; i++) {
for (int j=0 ; j<nt ; j++)
auxi.set(j) = 0 ;
auxi.set(i) = 1 ;
Tbl sortie (passage->inverse(auxi)) ;
for (int j=0 ; j<nt ; j++)
inv->set(j,i) = sortie(j) ;
}
}
// Fin du calcul des matrices...
for (int l=0 ; l<so.mg->get_nzone() ; l++)
if (so.c_cf->t[l]->get_etat() != ETATZERO)
for (int k=0 ; k<so.mg->get_np(l) ; k++)
for (int i=0 ; i<so.mg->get_nr(l) ; i++) {
Tbl coefs(nt) ;
coefs.set_etat_qcq() ;
for (int j=0 ; j<nt ; j++)
coefs.set(j) = (*so.c_cf)(l,k,j,i) ;
Tbl prod (nt) ;
prod.set_etat_qcq() ;
for (int j=0 ; j<nt ; j++)
prod.set(j) = 0 ;
if (sens) {
//calcul direct
for (int j=0 ; j<nt ; j++)
for (int jb=0 ; jb<nt ; jb++)
prod.set(j) += (*inv)(j, jb) * coefs(jb) ;
}
else {
//calcul inverse
for (int j=0 ; j<nt ; j++)
for (int jb=0 ; jb<nt ; jb++)
prod.set(j) += (*passage)(j, jb) * coefs(jb) ;
}
for (int j=0 ; j<nt ; j++)
so.c_cf->set(l,k,j,i) = prod(j) ;
}
// On met les bonnes bases :
int base_tet = so.base.get_base_t(0) ;
Base_val new_base (so.base) ;
if (sens) {
switch (base_tet) {
case T_COS_I:
new_base.set_base_t(T_CL_COS_I) ;
break ;
case T_SIN_I:
new_base.set_base_t(T_CL_SIN_I) ;
break ;
default:
cout << "Problem in rotate_propre_impair" << endl ;
abort() ;
break ;
}
}
else {
switch (base_tet) {
case T_CL_COS_I:
new_base.set_base_t(T_COS_I) ;
break ;
case T_CL_SIN_I:
new_base.set_base_t(T_SIN_I) ;
break ;
default:
cout << "Problem in rotate_propre_impair" << endl ;
abort() ;
break ;
}
}
so.c_cf->base = new_base ;
so.base = new_base ;
}
void Valeur::val_propre_1d() {
// On recupere la base en tet :
int nz = mg->get_nzone() ;
int base_tet = base.get_base_t(0) ;
// On verifie que c'est bien la meme partout...
for (int i=1 ; i<nz ; i++)
assert (base_tet == base.get_base_t(i)) ;
switch (base_tet) {
case T_COS_P:
rotate_propre_pair (*this, true) ;
break ;
case T_SIN_P:
rotate_propre_pair (*this, true) ;
break ;
case T_COS_I:
rotate_propre_impair (*this, true) ;
break ;
case T_SIN_I:
rotate_propre_impair (*this, true) ;
break ;
case T_CL_COS_P:
break ;
case T_CL_SIN_P:
break ;
case T_CL_COS_I:
break ;
case T_CL_SIN_I:
break ;
default:
cout << "Unknown basis in Valeur::val_propre_1d" << endl ;
abort() ;
break ;
}
}
void Valeur::val_propre_1d_i() {
// On recupere la base en tet :
int nz = mg->get_nzone() ;
int base_tet = base.get_base_t(0) ;
// On verifie que c'est bien la meme partout...
for (int i=1 ; i<nz ; i++)
assert (base_tet == base.get_base_t(i)) ;
switch (base_tet) {
case T_CL_COS_P:
rotate_propre_pair (*this, false) ;
break ;
case T_CL_SIN_P:
rotate_propre_pair (*this, false) ;
break ;
case T_CL_COS_I:
rotate_propre_impair (*this, false) ;
break ;
case T_CL_SIN_I:
rotate_propre_impair (*this, false) ;
break ;
case T_COS_P:
break ;
case T_SIN_P:
break ;
case T_COS_I:
break ;
case T_SIN_I:
break ;
default:
cout << "Unknown basis in Valeur::val_propre_1d_i" << endl ;
abort() ;
break ;
}
}
}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:322
Lorene::Base_val::set_base_t
void set_base_t(int base_t)
Sets the expansion basis for functions in all domains.
Definition base_val.C:195
Lorene::Base_val::get_base_t
int get_base_t(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
Definition base_val.h:411
Lorene::Matrice
Matrix handling.
Definition matrice.h:152
Lorene::Matrice::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition matrice.C:175
Lorene::Matrice::set
double & set(int j, int i)
Read/write of a particuliar element.
Definition matrice.h:277
Lorene::Matrice::inverse
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition matrice.C:424
Lorene::Matrice::set_band
void set_band(int up, int low) const
Calculate the band storage of *std.
Definition matrice.C:364
Lorene::Matrice::set_lu
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition matrice.C:392
Lorene::Matrice::vect_propre
Matrice vect_propre() const
Returns the eigenvectors of the matrix, calculated using LAPACK.
Definition matrice.C:507
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::Mtbl_cf::base
Base_val base
Bases of the spectral expansions.
Definition mtbl_cf.h:200
Lorene::Mtbl_cf::set
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:294
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::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::get_etat
int get_etat() const
Gives the logical state.
Definition tbl.h:394
Lorene::Tbl::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::set_etat_cf_qcq
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
Definition valeur.C:712
Lorene::Valeur::val_propre_1d_i
void val_propre_1d_i()
Inverse transformation of val_propre_1d.
Definition valeur_val_propre_1d.C:358
Lorene::Valeur::Valeur
Valeur(const Mg3d &mgrid)
Constructor.
Definition valeur.C:200
Lorene::Valeur::rotate_propre_pair
friend void rotate_propre_pair(Valeur &, bool)
Friend fonction.
Definition valeur_val_propre_1d.C:66
Lorene::Valeur::mg
const Mg3d * mg
Multi-grid Mgd3 on which this is defined.
Definition valeur.h:292
Lorene::Valeur::rotate_propre_impair
friend void rotate_propre_impair(Valeur &, bool)
Friend fonction.
Definition valeur_val_propre_1d.C:196
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::base
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
Lorene::Valeur::val_propre_1d
void val_propre_1d()
Set the basis to the eigenvalues of .
Definition valeur_val_propre_1d.C:322
T_CL_COS_I
#define T_CL_COS_I
CL of odd cosines.
Definition type_parite.h:232
T_CL_SIN_I
#define T_CL_SIN_I
CL of odd sines.
Definition type_parite.h:234
T_SIN_P
#define T_SIN_P
dev. sin seulement, harmoniques paires
Definition type_parite.h:202
T_CL_SIN_P
#define T_CL_SIN_P
CL of even sines.
Definition type_parite.h:230
T_COS_P
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition type_parite.h:200
T_CL_COS_P
#define T_CL_COS_P
CL of even cosines.
Definition type_parite.h:228
T_SIN_I
#define T_SIN_I
dev. sin seulement, harmoniques impaires
Definition type_parite.h:206
T_COS_I
#define T_COS_I
dev. cos seulement, harmoniques impaires
Definition type_parite.h:204
Lorene
Lorene prototypes.
Definition app_hor.h:64