scalar_raccord_externe.C
/*
* Copyright (c) 2003 Eric Gourgoulhon & Jerome Novak
*
* Copyright (c) 2001 Philippe Grandclement (for preceding Cmp version)
*
*
* 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 scalar_raccord_externe_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_raccord_externe.C,v 1.4 2014/10/13 08:53:47 j_novak Exp $" ;
/*
* $Id: scalar_raccord_externe.C,v 1.4 2014/10/13 08:53:47 j_novak Exp $
* $Log: scalar_raccord_externe.C,v $
* Revision 1.4 2014/10/13 08:53:47 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.3 2014/10/06 15:16:16 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.2 2003/09/25 09:22:33 j_novak
* Added a #include<math.h>
*
* Revision 1.1 2003/09/25 08:58:10 e_gourgoulhon
* First version.
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_raccord_externe.C,v 1.4 2014/10/13 08:53:47 j_novak Exp $
*
*/
//standard
#include <cstdlib>
#include <cmath>
// LORENE
#include "matrice.h"
#include "tensor.h"
#include "proto.h"
namespace Lorene {
int cnp (int n, int p) ;
// Fait le raccord dans la zec ...
// Suppose (pour le moment, le meme nbre de points sur les angles ...)
// et que la zone precedente est une coquille
void Scalar::raccord_externe(int power, int nbre, int lmax) {
va.coef() ;
va.ylm() ;
Base_val base_devel (va.base) ;
int base_r, m_quant, l_quant ;
// Confort :
int zone = mp->get_mg()->get_nzone()-2 ;
int nt = mp->get_mg()->get_nt(zone) ;
int np = mp->get_mg()->get_np(zone) ;
int nr = mp->get_mg()->get_nr(zone) ;
// Le mapping doit etre affine :
const Map_af* map = dynamic_cast<const Map_af*>(mp) ;
if (map == 0x0) {
cout << "Le mapping doit etre affine" << endl ;
abort() ;
}
// Mappinhg en r
double alpha = map->get_alpha()[zone] ;
double beta = map->get_beta()[zone] ;
// Mapping en 1/r
double new_alpha = -alpha/(beta*beta-alpha*alpha) ;
double new_beta = beta/(beta*beta-alpha*alpha) ;
// Mapping dans la zec :
double alpha_zec = map->get_alpha()[zone+1] ;
// Maintenant on construit les matrices de passage :
// Celle de ksi a T
Matrice tksi (nbre, nbre) ;
tksi.set_etat_qcq() ;
// Premier polynome
tksi.set(0, 0) = sqrt(double(2)) ;
for (int i=1 ; i<nbre ; i++)
tksi.set(0, i) = 0 ;
//Second polynome
tksi.set(1, 0) = 0 ;
tksi.set(1, 1) = sqrt(double(2)) ;
for (int i=2 ; i<nbre ; i++)
tksi.set(1, i) = 0 ;
// On recurre :
for (int lig=2 ; lig<nbre ; lig++) {
tksi.set(lig, 0) = -tksi(lig-2, 0) ;
for (int col=1 ; col<nbre ; col++)
tksi.set(lig, col) = 2*tksi(lig-1, col-1)-tksi(lig-2, col) ;
}
// Celle de u/new_alpha a ksi :
Matrice ksiu (nbre, nbre) ;
ksiu.set_etat_qcq() ;
for (int lig=0 ; lig<nbre ; lig++) {
for (int col=0 ; col<=lig ; col++)
ksiu.set(lig, col) = cnp(lig, col)*
pow(-new_beta/new_alpha, lig-col) ;
for (int col = lig+1 ; col<nbre ; col++)
ksiu.set(lig, col) = 0 ;
}
// La matrice totale :
Matrice tu (nbre, nbre) ;
tu.set_etat_qcq() ;
double somme ;
for (int lig=0 ; lig<nbre ; lig++)
for (int col=0 ; col<nbre ; col++) {
somme = 0 ;
for (int m=0 ; m<nbre ; m++)
somme += tksi(lig, m)*ksiu(m, col) ;
tu.set(lig, col) = somme ;
}
// On calcul les coefficients de u^n dans la zec
Tbl coef_u (nbre+lmax, nr) ;
coef_u.set_etat_qcq() ;
int* dege = new int [3] ;
dege[0] = 1 ; dege[1] = 1 ; dege[2] = nr ;
double* ti = new double [nr] ;
for (int puiss=0 ; puiss<nbre+lmax ; puiss++) {
for (int i=0 ; i<nr ; i++)
ti[i] = pow(-cos(M_PI*i/(nr-1))-1, puiss) ;
cfrcheb (dege, dege, ti, dege, ti) ;
for (int i=0 ; i<nr ; i++)
coef_u.set(puiss, i) = ti[i] ;
}
// Avant d entrer dans la boucle :
dege[2] = nbre ;
double *coloc = new double[nbre] ;
double *auxi = new double [1] ;
Tbl coef_zec (np+2, nt, nr) ;
coef_zec.annule_hard() ;
// Boucle sur les harmoniques :
for (int k=0 ; k<np+2 ; k++)
for (int j=0 ; j<nt ; j++)
if (nullite_plm (j, nt, k, np, base_devel)==1) {
donne_lm (zone+2, zone+1, j, k, base_devel, m_quant,
l_quant, base_r) ;
if (l_quant <= lmax) {
// On bosse :
// On recupere les valeus aux points de colocation en 1/r :
double ksi, air ;
for (int i=0 ; i<nbre ; i++) {
ksi = -cos(M_PI*i/(nbre-1)) ;
air = 1./(new_alpha*ksi+new_beta) ;
ksi = (air-beta)/alpha ;
for (int m=0 ; m<nr ; m++)
ti[m] = (*va.c_cf)(zone, k, j, m) ;
som_r_cheb (ti, nr, 1, 1, ksi, auxi) ;
coloc[i] = auxi[0]/
pow (-new_alpha*cos(M_PI*i/(nbre-1))+new_beta, power+l_quant);
}
cfrcheb (dege, dege, coloc, dege, coloc) ;
Tbl expansion (nbre) ;
expansion.set_etat_qcq() ;
for (int i=0 ; i<nbre ; i++) {
somme = 0 ;
for (int m=0 ; m<nbre ; m++)
somme += coloc[m]*tu(m, i) ;
expansion.set(i) = somme ;
}
for (int i=0 ; i<nr ; i++) {
somme = 0 ;
for (int m=0 ; m<nbre ; m++)
somme += coef_u(m+l_quant, i)*expansion(m)*
pow(alpha_zec, m+l_quant)/
pow(new_alpha, m) ;
coef_zec.set(k, j, i) = somme ;
}
}
}
va.set_etat_cf_qcq() ;
va.c_cf->set_etat_qcq() ;
va.c_cf->t[zone+1]->set_etat_qcq() ;
for (int k=0 ; k<np+2 ; k++)
for (int j=0 ; j<nt ; j++)
for (int i=0 ; i<nr ; i++)
va.c_cf->set(zone+1, k, j, i) = coef_zec(k, j, i) ;
set_dzpuis(power) ;
va.ylm_i() ;
delete[] auxi ;
delete [] dege ;
delete [] ti ;
delete [] coloc ;
}
}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:322
Lorene::Map_af
Affine radial mapping.
Definition map.h:2027
Lorene::Map_af::get_beta
const double * get_beta() const
Returns the pointer on the array beta.
Definition map_af.C:481
Lorene::Map_af::get_alpha
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition map_af.C:477
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::Scalar::cos
friend Scalar cos(const Scalar &)
Cosine.
Definition scalar_math.C:104
Lorene::Scalar::raccord_externe
void raccord_externe(int puis, int nbre, int lmax)
Matching of the external domain with the outermost shell.
Definition scalar_raccord_externe.C:66
Lorene::Scalar::set_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
Lorene::Scalar::va
Valeur va
The numerical value of the Scalar.
Definition scalar.h:405
Lorene::Scalar::sqrt
friend Scalar sqrt(const Scalar &)
Square root.
Definition scalar_math.C:263
Lorene::Scalar::pow
friend Scalar pow(const Scalar &, int)
Power .
Definition scalar_math.C:451
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:372
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::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Lorene
Lorene prototypes.
Definition app_hor.h:64