sym_tensor_trans_pde.C
/*
* Functions to solve various PDEs for a divergence-free symmetric tensor.
*
* (see file sym_tensor.h for documentation).
*
*/
/*
* Copyright (c) 2005-2006 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 version 2
* as published by the Free Software Foundation.
*
* 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 sym_tensor_trans_pde_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_trans_pde.C,v 1.16 2014/10/13 08:53:43 j_novak Exp $" ;
/*
* $Id: sym_tensor_trans_pde.C,v 1.16 2014/10/13 08:53:43 j_novak Exp $
* $Log: sym_tensor_trans_pde.C,v $
* Revision 1.16 2014/10/13 08:53:43 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.15 2014/10/06 15:13:19 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.14 2010/10/11 10:38:34 j_novak
* *** empty log message ***
*
* Revision 1.13 2010/10/11 10:23:03 j_novak
* Removed methods Sym_tensor_trans::solve_hrr() and Sym_tensor_trans::set_WX_det_one(), as they are no longer relevant.
*
* Revision 1.12 2006/09/05 15:38:45 j_novak
* The fuctions sol_Dirac... are in a seperate file, with new parameters to
* control the boundary conditions.
*
* Revision 1.11 2006/08/31 12:13:22 j_novak
* Added an argument of type Param to Sym_tensor_trans::sol_Dirac_A().
*
* Revision 1.10 2006/06/28 07:48:26 j_novak
* Better treatment of some null cases.
*
* Revision 1.9 2006/06/21 15:42:47 j_novak
* Minor changes.
*
* Revision 1.8 2006/06/20 12:07:15 j_novak
* Improved execution speed for sol_Dirac_tildeB...
*
* Revision 1.7 2006/06/14 10:04:21 j_novak
* New methods sol_Dirac_l01, set_AtB_det_one and set_AtB_trace_zero.
*
* Revision 1.6 2006/06/13 13:30:12 j_novak
* New members sol_Dirac_A and sol_Dirac_tildeB (see documentation).
*
* Revision 1.5 2006/06/12 13:37:23 j_novak
* Added bounds in l (multipolar momentum) for Sym_tensor_trans::solve_hrr.
*
* Revision 1.4 2005/11/28 14:45:17 j_novak
* Improved solution of the Poisson tensor equation in the case of a transverse
* tensor.
*
* Revision 1.3 2005/11/24 14:07:54 j_novak
* Use of Matrice::annule_hard()
*
* Revision 1.2 2005/11/24 09:24:25 j_novak
* Corrected some missing references.
*
* Revision 1.1 2005/09/16 13:58:11 j_novak
* New Poisson solver for a Sym_tensor_trans.
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_trans_pde.C,v 1.16 2014/10/13 08:53:43 j_novak Exp $
*
*/
// C headers
#include <cassert>
#include <cmath>
// Lorene headers
#include "tensor.h"
#include "diff.h"
#include "proto.h"
#include "param.h"
namespace Lorene {
Sym_tensor_trans Sym_tensor_trans::poisson(const Scalar* h_guess) const {
// All this has a meaning only for spherical components...
assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
//## ... and affine mapping, for the moment!
const Map_af* mpaff = dynamic_cast<const Map_af*>(mp) ;
assert( mpaff!= 0x0) ;
Sym_tensor_trans resu(*mp, *triad, *met_div) ;
const Mg3d& gri = *mp->get_mg() ;
int np = gri.get_np(0) ;
int nt = gri.get_nt(0) ;
assert (nt > 4) ;
if (np == 1) {
int nz = gri.get_nzone() ;
double* bornes = new double[nz+1] ;
const double* alp = mpaff->get_alpha() ;
const double* bet = mpaff->get_beta() ;
for (int lz=0; lz<nz; lz++) {
assert (gri.get_np(lz) == np) ;
assert (gri.get_nt(lz) == nt) ;
switch (gri.get_type_r(lz)) {
case RARE: {
bornes[lz] = bet[lz] ;
break ;
}
case FIN: {
bornes[lz] = bet[lz] - alp[lz] ;
break ;
}
case UNSURR: {
bornes[lz] = double(1) / ( bet[lz] - alp[lz] ) ;
break ;
}
default: {
cout << "Sym_tensor_trans::poisson() : problem with the grid!"
<< endl ;
abort() ;
break ;
}
}
}
if (gri.get_type_r(nz-1) == UNSURR)
bornes[nz] = 1./(alp[nz-1] + bet[nz-1]) ;
else
bornes[nz] = alp[nz-1] + bet[nz-1] ;
const Mg3d& gr2 = *gri.get_non_axi() ;
Map_af mp2(gr2, bornes) ;
int np2 = ( np > 3 ? np : 4 ) ;
Sym_tensor sou_cart(mp2, CON, mp2.get_bvect_spher()) ;
for (int l=1; l<=3; l++)
for (int c=l; c<=3; c++) {
switch (this->operator()(l,c).get_etat() ) {
case ETATZERO: {
sou_cart.set(l,c).set_etat_zero() ;
break ;
}
case ETATUN: {
sou_cart.set(l,c).set_etat_one() ;
break ;
}
case ETATQCQ : {
sou_cart.set(l,c).allocate_all() ;
for (int lz=0; lz<nz; lz++)
for (int k=0; k<np2; k++)
for (int j=0; j<nt; j++)
for(int i=0; i<gr2.get_nr(lz); i++)
sou_cart.set(l,c).set_grid_point(lz, k, j, i)
= this->operator()(l,c).val_grid_point(lz, 0, j, i) ;
break ;
}
default: {
cout <<
"Sym_tensor_trans::poisson() : source in undefined state!"
<< endl ;
abort() ;
break ;
}
}
sou_cart.set(l,c).set_dzpuis(this->operator()(l,c).get_dzpuis()) ;
}
sou_cart.std_spectral_base() ;
sou_cart.change_triad(mp2.get_bvect_cart()) ;
Sym_tensor res_cart(mp2, CON, mp2.get_bvect_cart()) ;
for (int i=1; i<=3; i++)
for(int j=i; j<=3; j++)
res_cart.set(i,j) = sou_cart(i,j).poisson() ;
res_cart.change_triad(mp2.get_bvect_spher()) ;
Scalar res_A(*mp) ; Scalar big_A = res_cart.compute_A() ;
Scalar res_B(*mp) ; Scalar big_B = res_cart.compute_tilde_B_tt() ;
switch (big_A.get_etat() ) {
case ETATZERO: {
res_A.set_etat_zero() ;
break ;
}
case ETATUN : {
res_A.set_etat_one() ;
break ;
}
case ETATQCQ : {
res_A.allocate_all() ;
for (int lz=0; lz<nz; lz++)
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++)
for(int i=0; i<gri.get_nr(lz); i++)
res_A.set_grid_point(lz, k, j, i)
= big_A.val_grid_point(lz, k, j, i) ;
break ;
}
default: {
cout <<
"Sym_tensor_trans::poisson() : res_A in undefined state!"
<< endl ;
abort() ;
break ;
}
}
res_A.set_spectral_base(big_A.get_spectral_base()) ;
int dzA = big_A.get_dzpuis() ;
res_A.set_dzpuis(dzA) ;
switch (big_B.get_etat() ) {
case ETATZERO: {
res_B.set_etat_zero() ;
break ;
}
case ETATUN : {
res_B.set_etat_one() ;
break ;
}
case ETATQCQ : {
res_B.allocate_all() ;
for (int lz=0; lz<nz; lz++)
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++)
for(int i=0; i<gri.get_nr(lz); i++)
res_B.set_grid_point(lz, k, j, i)
= big_B.val_grid_point(lz, k, j, i) ;
break ;
}
default: {
cout <<
"Sym_tensor_trans::poisson() : res_B in undefined state!"
<< endl ;
abort() ;
break ;
}
}
res_B.set_spectral_base(big_B.get_spectral_base()) ;
int dzB = big_B.get_dzpuis() ;
res_B.set_dzpuis(dzB) ;
resu.set_AtBtt_det_one(res_A, res_B, h_guess) ;
delete [] bornes ;
}
else {
assert (np >=4) ;
Sym_tensor_trans sou_cart = *this ;
sou_cart.change_triad(mp->get_bvect_cart()) ;
Sym_tensor res_cart(*mp, CON, mp->get_bvect_cart()) ;
for (int i=1; i<=3; i++)
for(int j=i; j<=3; j++)
res_cart.set(i,j) = sou_cart(i,j).poisson() ;
res_cart.change_triad(*triad) ;
resu.set_AtBtt_det_one(res_cart.compute_A(), res_cart.compute_tilde_B_tt(), h_guess) ;
}
#ifndef NDEBUG
Vector dive = resu.divergence(*met_div) ;
dive.dec_dzpuis(2) ;
maxabs(dive, "Sym_tensor_trans::poisson : divergence of the solution") ;
#endif
return resu ;
}
}
Lorene::Base_vect_spher
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
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::Mg3d
Multi-domain grid.
Definition grilles.h:273
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::get_type_r
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:474
Lorene::Mg3d::get_non_axi
const Mg3d * get_non_axi() const
Returns the pointer on the grid which has at least 4 points in the direction and at least 5 in the ...
Definition mg3d.C:615
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::set_etat_one
void set_etat_one()
Sets the logical state to ETATUN (one).
Definition scalar.C:334
Lorene::Scalar::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition scalar.h:557
Lorene::Scalar::val_grid_point
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition scalar.h:637
Lorene::Scalar::set_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
Lorene::Scalar::allocate_all
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition scalar.C:365
Lorene::Scalar::get_spectral_base
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Definition scalar.h:1294
Lorene::Scalar::set_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
Lorene::Scalar::set_grid_point
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition scalar.h:684
Lorene::Scalar::set_spectral_base
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va.
Definition scalar.C:797
Lorene::Sym_tensor_trans::set_AtBtt_det_one
void set_AtBtt_det_one(const Scalar &a_in, const Scalar &tbtt_in, const Scalar *h_prev=0x0, Param *par_bc=0x0, Param *par_mat=0x0, double precis=1.e-14, int it_max=100)
Assigns the derived member A and computes from its TT-part (see Sym_tensor::compute_tilde_B_tt() ).
Definition sym_tensor_trans_aux.C:137
Lorene::Sym_tensor_trans::Sym_tensor_trans
Sym_tensor_trans(const Map &map, const Base_vect &triad_i, const Metric &met)
Standard constructor.
Definition sym_tensor_trans.C:114
Lorene::Sym_tensor_trans::poisson
Sym_tensor_trans poisson(const Scalar *h_guess=0x0) const
Computes the solution of a tensorial transverse Poisson equation with *this as a source:
Definition sym_tensor_trans_pde.C:99
Lorene::Sym_tensor_trans::met_div
const Metric *const met_div
Metric with respect to which the divergence and the trace are defined.
Definition sym_tensor.h:614
Lorene::Sym_tensor::compute_A
const Scalar & compute_A(bool output_ylm=true, Param *par=0x0) const
Gives the field A (see member p_aaa ).
Definition sym_tensor_aux.C:316
Lorene::Sym_tensor::compute_tilde_B_tt
Scalar compute_tilde_B_tt(bool output_ylm=true, Param *par=0x0) const
Gives the field (see member p_tilde_b ) associated with the TT-part of the Sym_tensor .
Definition sym_tensor_aux.C:478
Lorene::Sym_tensor::Sym_tensor
Sym_tensor(const Map &map, const Itbl &tipe, const Base_vect &triad_i)
Standard constructor.
Definition sym_tensor.C:142
Lorene::Sym_tensor::divergence
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition sym_tensor.C:349
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene::Tensor::dec_dzpuis
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:808
Lorene::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Lorene::Tensor::change_triad
virtual void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition tensor_change_triad.C:77
Lorene::Tensor::operator()
const Scalar & operator()(const Itbl &ind) const
Returns the value of a component (read-only version).
Definition tensor.C:798
Lorene::Tensor::set
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
Lorene::Tensor::triad
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tensor.h:303
Lorene::Tensor::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition tensor.C:926
Lorene::maxabs
Tbl maxabs(const Tensor &aa, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maxima in each domain of the absolute values of the tensor components.
Definition tensor_calculus_ext.C:606
Lorene
Lorene prototypes.
Definition app_hor.h:64