tenseur_pde.C
/*
* Copyright (c) 2000-2001 Eric Gourgoulhon
* Copyright (c) 2000-2001 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 tenseur_pde_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_pde.C,v 1.7 2014/10/13 08:53:42 j_novak Exp $" ;
/*
* $Id: tenseur_pde.C,v 1.7 2014/10/13 08:53:42 j_novak Exp $
* $Log: tenseur_pde.C,v $
* Revision 1.7 2014/10/13 08:53:42 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.6 2006/06/01 12:47:54 p_grandclement
* update of the Bin_ns_bh project
*
* Revision 1.5 2005/08/30 08:35:13 p_grandclement
* Addition of the Tau version of the vectorial Poisson equation for the Tensors
*
* Revision 1.4 2003/10/03 15:58:51 j_novak
* Cleaning of some headers
*
* Revision 1.3 2002/08/07 16:14:11 j_novak
* class Tenseur can now also handle tensor densities, this should be transparent to older codes
*
* Revision 1.2 2002/07/09 16:46:23 p_grandclement
* The Param in the case of an affine mapping is now 0x0 and not deleted
* (I wonder why it was working before)
*
* Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
* LORENE
*
* Revision 2.13 2000/10/04 14:58:32 eric
* Ajout de shift.set_etat_qcq() avant l'affection de shift.
*
* Revision 2.12 2000/09/27 15:07:22 eric
* Traitement source nulle dans poisson_vect.
*
* Revision 2.11 2000/05/22 15:48:18 phil
* modification de Oohara pour passer avec dzpuis == 2
* pardon 3
*
* Revision 2.10 2000/05/22 15:00:46 phil
* Modification de la methode de Shibata :
* on doit passer une source en r^4 et l equation scalaire est alors resolue en utilisant l'algotihme avec dzpuis == 3
*
* Revision 2.9 2000/03/10 15:51:42 eric
* Traitement dzpuis de source_scal.
*
* Revision 2.8 2000/03/08 10:21:05 eric
* Appel de delete sur le Param* retourne par Map_et::donne_para_poisson_vect[
* lorsqu'il n'est plus utilise (correction Memory leak).
*
* Revision 2.7 2000/03/07 16:53:42 eric
* *** empty log message ***
*
* Revision 2.6 2000/03/07 15:43:32 phil
* gestion des cas dzpuis ==4
*
* Revision 2.5 2000/02/21 12:55:09 eric
* Traitement des triades.
*
* Revision 2.4 2000/02/16 17:13:05 eric
* Correction
* mp->donne_para_poisson_vect(para, 4)
* devient
* mp->donne_para_poisson_vect(para, 3)
* dans Tenseur::poisson_vect
* /
*
* Revision 2.3 2000/02/15 10:26:49 phil
* le calcul de sol n'appelle plus Map::poisson_vect mais est fait dans
* Tenseur::poisson_vect (respectivement poisson_vect_oohara)
*
* Revision 2.2 2000/02/09 19:32:58 eric
* La triade de decomposition est desormais passee en argument des constructeurs.
*
* Revision 2.1 2000/02/09 10:01:32 phil
* ajout version oohara
*
* Revision 2.0 2000/01/21 12:58:57 phil
* *** empty log message ***
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_pde.C,v 1.7 2014/10/13 08:53:42 j_novak Exp $
*
*/
// Header Lorene:
#include "param.h"
#include "tenseur.h"
//-----------------------------------//
// Vectorial Poisson equation //
//-----------------------------------//
// Version avec parametres
// -----------------------
namespace Lorene {
void Tenseur::poisson_vect(double lambda, Param& para, Tenseur& shift
, Tenseur& vecteur, Tenseur& scalaire) const {
assert (lambda != -1) ;
// Verifications d'usage ...
assert (valence == 1) ;
assert (shift.get_valence() == 1) ;
assert (shift.get_type_indice(0) == type_indice(0)) ;
assert (vecteur.get_valence() == 1) ;
assert (vecteur.get_type_indice(0) == type_indice(0)) ;
assert (scalaire.get_valence() == 0) ;
assert (etat != ETATNONDEF) ;
// Nothing to do if the source is zero
if (etat == ETATZERO) {
shift.set_etat_zero() ;
vecteur.set_etat_qcq() ;
for (int i=0; i<3; i++) {
vecteur.set(i) = 0 ;
}
scalaire.set_etat_qcq() ;
scalaire.set() = 0 ;
return ;
}
for (int i=0 ; i<3 ; i++)
assert ((*this)(i).check_dzpuis(4)) ;
// On construit le tableau contenant le terme P_i ...
for (int i=0 ; i<3 ; i++) {
Param* par = mp->donne_para_poisson_vect(para, i) ;
(*this)(i).poisson(*par, vecteur.set(i)) ;
if (par != 0x0)
delete par ;
}
vecteur.set_triad( *triad ) ;
// Equation de Poisson scalaire :
Tenseur source_scal (-skxk(*this)) ;
assert (source_scal().check_dzpuis(3)) ;
Param* par = mp->donne_para_poisson_vect(para, 3) ;
source_scal().poisson(*par, scalaire.set()) ;
if (par !=0x0)
delete par ;
// On construit le tableau contenant le terme d xsi / d x_i ...
Tenseur auxiliaire(scalaire) ;
Tenseur dxsi (auxiliaire.gradient()) ;
dxsi.dec2_dzpuis() ;
// On construit le tableau contenant le terme x_k d P_k / d x_i
Tenseur dp (skxk(vecteur.gradient())) ;
dp.dec_dzpuis() ;
// Il ne reste plus qu'a tout ranger dans P :
// The final computation is done component by component because
// d_khi and x_d_w are covariant comp. whereas w_shift is
// contravariant
shift.set_etat_qcq() ;
for (int i=0 ; i<3 ; i++)
shift.set(i) = (lambda+2)/2/(lambda+1) * vecteur(i)
- (lambda/2/(lambda+1)) * (dxsi(i) + dp(i)) ;
shift.set_triad( *(vecteur.triad) ) ;
}
// Version sans parametres
// -----------------------
Tenseur Tenseur::poisson_vect(double lambda, Tenseur& vecteur,
Tenseur& scalaire) const {
Param bidon ;
Tenseur resu(*mp, valence, type_indice, triad, metric, poids) ;
resu.set_etat_qcq() ;
poisson_vect(lambda, bidon, resu, vecteur, scalaire) ;
return resu ;
}
//-----------------------------------//
// Vectorial Poisson equation //
// using Oohara scheme //
//-----------------------------------//
// Version avec parametres
// -----------------------
void Tenseur::poisson_vect_oohara(double lambda, Param& para, Tenseur& shift,
Tenseur& chi) const {
// Ne marche pas pour lambda =-1
assert (lambda != -1) ;
// Verifications d'usage ...
assert (valence == 1) ;
assert (shift.get_valence() == 1) ;
assert (shift.get_type_indice(0) == type_indice(0)) ;
assert (chi.get_valence() == 0) ;
assert (etat != ETATNONDEF) ;
// Nothing to do if the source is zero
if (etat == ETATZERO) {
shift.set_etat_zero() ;
chi.set_etat_qcq() ;
chi.set() = 0 ;
return ;
}
for (int i=0 ; i<3 ; i++)
assert ((*this)(i).check_dzpuis(3) ||
(*this)(i).check_dzpuis(4)) ;
Tenseur copie(*this) ;
copie.dec2_dzpuis() ;
if ((*this)(0).check_dzpuis(4))
copie.dec2_dzpuis() ;
else
copie.dec_dzpuis() ;
Tenseur source_scal(contract(copie.gradient(), 0, 1)/(1.+lambda)) ;
source_scal.inc2_dzpuis() ;
Param* par = mp->donne_para_poisson_vect(para, 3) ;
source_scal().poisson(*par, chi.set());
if (par !=0x0)
delete par ;
Tenseur source_vect(*this) ;
if ((*this)(0).check_dzpuis(4))
source_vect.dec_dzpuis() ;
Tenseur chi_grad (chi.gradient()) ;
chi_grad.inc_dzpuis() ;
for (int i=0 ; i<3 ; i++)
source_vect.set(i) -= lambda*chi_grad(i) ;
assert( *(source_vect.triad) == *((chi.gradient()).get_triad()) ) ;
if (shift.get_etat() == ETATZERO) {
shift.set_etat_qcq() ;
for (int i=0 ; i<3 ; i++)
shift.set(i) = 0 ;
}
for (int i=0 ; i<3 ; i++) {
par = mp->donne_para_poisson_vect(para, i) ;
source_vect(i).poisson(*par, shift.set(i)) ;
if (par !=0x0)
delete par ;
}
shift.set_triad( *(source_vect.triad) ) ;
}
// Version sans parametres
// -----------------------
Tenseur Tenseur::poisson_vect_oohara(double lambda, Tenseur& scalaire) const {
Param bidon ;
Tenseur resu(*mp, valence, type_indice, triad, metric, poids) ;
resu.set_etat_qcq() ;
poisson_vect_oohara(lambda, bidon, resu, scalaire) ;
return resu ;
}
//---------------------------------------------//
// Vectorial Poisson equation TAU method//
//---------------------------------------------//
// Version avec parametres
// -----------------------
void Tenseur::poisson_vect_tau(double lambda, Param& para, Tenseur& shift
, Tenseur& vecteur, Tenseur& scalaire) const {
assert (lambda != -1) ;
// Verifications d'usage ...
assert (valence == 1) ;
assert (shift.get_valence() == 1) ;
assert (shift.get_type_indice(0) == type_indice(0)) ;
assert (vecteur.get_valence() == 1) ;
assert (vecteur.get_type_indice(0) == type_indice(0)) ;
assert (scalaire.get_valence() == 0) ;
assert (etat != ETATNONDEF) ;
// Nothing to do if the source is zero
if (etat == ETATZERO) {
shift.set_etat_zero() ;
vecteur.set_etat_qcq() ;
for (int i=0; i<3; i++) {
vecteur.set(i) = 0 ;
}
scalaire.set_etat_qcq() ;
scalaire.set() = 0 ;
return ;
}
for (int i=0 ; i<3 ; i++)
assert ((*this)(i).check_dzpuis(4)) ;
// On construit le tableau contenant le terme P_i ...
for (int i=0 ; i<3 ; i++) {
Param* par = mp->donne_para_poisson_vect(para, i) ;
(*this)(i).poisson_tau(*par, vecteur.set(i)) ;
if (par != 0x0)
delete par ;
}
vecteur.set_triad( *triad ) ;
// Equation de Poisson scalaire :
Tenseur source_scal (-skxk(*this)) ;
assert (source_scal().check_dzpuis(3)) ;
Param* par = mp->donne_para_poisson_vect(para, 3) ;
source_scal().poisson_tau(*par, scalaire.set()) ;
if (par !=0x0)
delete par ;
// On construit le tableau contenant le terme d xsi / d x_i ...
Tenseur auxiliaire(scalaire) ;
Tenseur dxsi (auxiliaire.gradient()) ;
dxsi.dec2_dzpuis() ;
// On construit le tableau contenant le terme x_k d P_k / d x_i
Tenseur dp (skxk(vecteur.gradient())) ;
dp.dec_dzpuis() ;
// Il ne reste plus qu'a tout ranger dans P :
// The final computation is done component by component because
// d_khi and x_d_w are covariant comp. whereas w_shift is
// contravariant
shift.set_etat_qcq() ;
for (int i=0 ; i<3 ; i++)
shift.set(i) = (lambda+2)/2/(lambda+1) * vecteur(i)
- (lambda/2/(lambda+1)) * (dxsi(i) + dp(i)) ;
shift.set_triad( *(vecteur.triad) ) ;
}
// Version sans parametres
// -----------------------
Tenseur Tenseur::poisson_vect_tau(double lambda, Tenseur& vecteur,
Tenseur& scalaire) const {
Param bidon ;
Tenseur resu(*mp, valence, type_indice, triad, metric, poids) ;
resu.set_etat_qcq() ;
poisson_vect_tau(lambda, bidon, resu, vecteur, scalaire) ;
return resu ;
}
//-----------------------------------//
// Vectorial Poisson equation //
// using Oohara scheme //
//-----------------------------------//
// Version avec parametres
// -----------------------
void Tenseur::poisson_vect_oohara_tau(double lambda, Param& para, Tenseur& shift,
Tenseur& chi) const {
// Ne marche pas pour lambda =-1
assert (lambda != -1) ;
// Verifications d'usage ...
assert (valence == 1) ;
assert (shift.get_valence() == 1) ;
assert (shift.get_type_indice(0) == type_indice(0)) ;
assert (chi.get_valence() == 0) ;
assert (etat != ETATNONDEF) ;
// Nothing to do if the source is zero
if (etat == ETATZERO) {
shift.set_etat_zero() ;
chi.set_etat_qcq() ;
chi.set() = 0 ;
return ;
}
for (int i=0 ; i<3 ; i++)
assert ((*this)(i).check_dzpuis(3) ||
(*this)(i).check_dzpuis(4)) ;
Tenseur copie(*this) ;
copie.dec2_dzpuis() ;
if ((*this)(0).check_dzpuis(4))
copie.dec2_dzpuis() ;
else
copie.dec_dzpuis() ;
Tenseur source_scal(contract(copie.gradient(), 0, 1)/(1.+lambda)) ;
source_scal.inc2_dzpuis() ;
Param* par = mp->donne_para_poisson_vect(para, 3) ;
source_scal().poisson_tau(*par, chi.set());
if (par !=0x0)
delete par ;
Tenseur source_vect(*this) ;
if ((*this)(0).check_dzpuis(4))
source_vect.dec_dzpuis() ;
Tenseur chi_grad (chi.gradient()) ;
chi_grad.inc_dzpuis() ;
for (int i=0 ; i<3 ; i++)
source_vect.set(i) -= lambda*chi_grad(i) ;
assert( *(source_vect.triad) == *((chi.gradient()).get_triad()) ) ;
for (int i=0 ; i<3 ; i++) {
par = mp->donne_para_poisson_vect(para, i) ;
source_vect(i).poisson_tau(*par, shift.set(i)) ;
if (par !=0x0)
delete par ;
}
shift.set_triad( *(source_vect.triad) ) ;
}
// Version sans parametres
// -----------------------
Tenseur Tenseur::poisson_vect_oohara_tau(double lambda, Tenseur& scalaire) const {
Param bidon ;
Tenseur resu(*mp, valence, type_indice, triad, metric, poids) ;
resu.set_etat_qcq() ;
poisson_vect_oohara_tau(lambda, bidon, resu, scalaire) ;
return resu ;
}
}
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Lorene::Tenseur::set
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
Lorene::Tenseur::get_type_indice
int get_type_indice(int i) const
Returns the type of the index number i .
Definition tenseur.h:726
Lorene::Tenseur::metric
const Metrique * metric
For tensor densities: the metric defining the conformal factor.
Definition tenseur.h:325
Lorene::Tenseur::mp
const Map *const mp
Reference mapping.
Definition tenseur.h:306
Lorene::Tenseur::inc_dzpuis
void inc_dzpuis()
dzpuis += 1 ;
Definition tenseur.C:1117
Lorene::Tenseur::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
Lorene::Tenseur::poisson_vect_oohara
void poisson_vect_oohara(double lambda, Param &par, Tenseur &shift, Tenseur &scal) const
Solves the vectorial Poisson equation .
Definition tenseur_pde.C:218
Lorene::Tenseur::Tenseur
Tenseur(const Map &map, const Metrique *met=0x0, double weight=0)
Constructor for a scalar field.
Definition tenseur.C:209
Lorene::Tenseur::dec2_dzpuis
void dec2_dzpuis()
dzpuis -= 2 ;
Definition tenseur.C:1130
Lorene::Tenseur::dec_dzpuis
void dec_dzpuis()
dzpuis -= 1 ;
Definition tenseur.C:1104
Lorene::Tenseur::triad
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tenseur.h:312
Lorene::Tenseur::gradient
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates).
Definition tenseur.C:1542
Lorene::Tenseur::contract
friend Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition tenseur_operateur.C:279
Lorene::Tenseur::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition tenseur.C:645
Lorene::Tenseur::type_indice
Itbl type_indice
Array of size valence contening the type of each index, COV for a covariant one and CON for a contrav...
Definition tenseur.h:318
Lorene::Tenseur::valence
int valence
Valence.
Definition tenseur.h:307
Lorene::Tenseur::inc2_dzpuis
void inc2_dzpuis()
dzpuis += 2 ;
Definition tenseur.C:1143
Lorene::Tenseur::etat
int etat
Logical state ETATZERO , ETATQCQ or ETATNONDEF.
Definition tenseur.h:321
Lorene::Tenseur::poids
double poids
For tensor densities: the weight.
Definition tenseur.h:323
Lorene::Tenseur::skxk
friend Tenseur skxk(const Tenseur &)
Contraction of the last index of (*this) with or , depending on the type of S .
Definition tenseur_operateur.C:580
Lorene::Tenseur::get_valence
int get_valence() const
Returns the valence.
Definition tenseur.h:710
Lorene::Tenseur::set_triad
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition tenseur.C:674
Lorene::Tenseur::poisson_vect
void poisson_vect(double lambda, Param &par, Tenseur &shift, Tenseur &vect, Tenseur &scal) const
Solves the vectorial Poisson equation : .
Definition tenseur_pde.C:118
Lorene::Tenseur::get_etat
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::poisson_tau
virtual void poisson_tau(const Cmp &source, Param &par, Cmp &uu) const =0
Computes the solution of a scalar Poisson equationwith a Tau method.
Lorene::poisson
virtual void poisson(const Cmp &source, Param &par, Cmp &uu) const =0
Computes the solution of a scalar Poisson equation.