LORENE
tenseur_pde_falloff.C
1 /*
2  * Methods of the class tenseur for solving vectorial Poisson equations
3  * with a falloff condition at the outer boundary
4  *
5  * (see file tenseur.h for documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 2004 Joshua A. Faber
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License version 2
16  * as published by the Free Software Foundation.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 char tenseur_pde_falloff_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_pde_falloff.C,v 1.3 2014/10/13 08:53:42 j_novak Exp $" ;
30 
31 /*
32  * $Id: tenseur_pde_falloff.C,v 1.3 2014/10/13 08:53:42 j_novak Exp $
33  * $Log: tenseur_pde_falloff.C,v $
34  * Revision 1.3 2014/10/13 08:53:42 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.2 2004/12/22 18:25:59 k_taniguchi
38  * Cahnge an argument of poisson_vect_falloff.
39  *
40  * Revision 1.1 2004/11/30 20:55:33 k_taniguchi
41  * *** empty log message ***
42  *
43  *
44  * $Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_pde_falloff.C,v 1.3 2014/10/13 08:53:42 j_novak Exp $
45  *
46  */
47 
48 // Header Lorene:
49 #include "map.h"
50 #include "cmp.h"
51 #include "param.h"
52 #include "tenseur.h"
53 
54  //-----------------------------------//
55  // Vectorial Poisson equation //
56  //-----------------------------------//
57 
58 // Version avec parametres
59 // -----------------------
60 namespace Lorene {
61 void Tenseur::poisson_vect_falloff(double lambda, Param& para, Tenseur& shift,
62  Tenseur& vecteur, Tenseur& scalaire,
63  int* k_falloff) const {
64  assert (lambda != -1) ;
65 
66  // Verifications d'usage ...
67  assert (valence == 1) ;
68  assert (shift.get_valence() == 1) ;
69  assert (shift.get_type_indice(0) == type_indice(0)) ;
70  assert (vecteur.get_valence() == 1) ;
71  assert (vecteur.get_type_indice(0) == type_indice(0)) ;
72  assert (scalaire.get_valence() == 0) ;
73  assert (etat != ETATNONDEF) ;
74 
75  // define int k_falloff[4] ;
76  // k_falloff[0,1,2] = falloff for vecteur_x,y,z
77  // k_falloff[3] = falloff for scalaire
78 
79  // Nothing to do if the source is zero
80  if (etat == ETATZERO) {
81 
82  shift.set_etat_zero() ;
83 
84  vecteur.set_etat_qcq() ;
85  for (int i=0; i<3; i++) {
86  vecteur.set(i) = 0 ;
87  }
88 
89  scalaire.set_etat_qcq() ;
90  scalaire.set() = 0 ;
91 
92  return ;
93  }
94 
95  // On construit le tableau contenant le terme P_i ...
96  for (int i=0 ; i<3 ; i++) {
97  Param* par = mp->donne_para_poisson_vect(para, i) ;
98 
99  (*this)(i).poisson_falloff(*par, vecteur.set(i),k_falloff[i]) ;
100 
101  if (par != 0x0)
102  delete par ;
103  }
104  vecteur.set_triad( *triad ) ;
105 
106  // Equation de Poisson scalaire :
107  Tenseur source_scal (-skxk(*this)) ;
108 
109  Param* par = mp->donne_para_poisson_vect(para, 3) ;
110 
111  source_scal().poisson_falloff(*par, scalaire.set(), k_falloff[3]) ;
112 
113  if (par !=0x0)
114  delete par ;
115 
116  // On construit le tableau contenant le terme d xsi / d x_i ...
117  Tenseur auxiliaire(scalaire) ;
118  Tenseur dxsi (auxiliaire.gradient()) ;
119 
120  // On construit le tableau contenant le terme x_k d P_k / d x_i
121  Tenseur dp (skxk(vecteur.gradient())) ;
122 
123  // Il ne reste plus qu'a tout ranger dans P :
124  // The final computation is done component by component because
125  // d_khi and x_d_w are covariant comp. whereas w_shift is
126  // contravariant
127 
128  shift.set_etat_qcq() ;
129 
130  for (int i=0 ; i<3 ; i++)
131  shift.set(i) = (lambda+2)/2/(lambda+1) * vecteur(i)
132  - (lambda/2/(lambda+1)) * (dxsi(i) + dp(i)) ;
133 
134  shift.set_triad( *(vecteur.triad) ) ;
135 
136 }
137 
138 
139 // Version sans parametres
140 // -----------------------
141 Tenseur Tenseur::poisson_vect_falloff(double lambda, Tenseur& vecteur,
142  Tenseur& scalaire, int* k_falloff) const {
143 
144  Param bidon ;
146  resu.set_etat_qcq() ;
147  poisson_vect_falloff(lambda, bidon, resu, vecteur, scalaire, k_falloff) ;
148  return resu ;
149 }
150 
151 }
Lorene::Tenseur::etat
int etat
Logical state ETATZERO , ETATQCQ or ETATNONDEF.
Definition: tenseur.h:321
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Tenseur::valence
int valence
Valence.
Definition: tenseur.h:307
Lorene::Tenseur::poids
double poids
For tensor densities: the weight.
Definition: tenseur.h:323
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::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::Tenseur
Tenseur(const Map &map, const Metrique *met=0x0, double weight=0)
Constructor for a scalar field.
Definition: tenseur.C:209
Lorene::Map::donne_para_poisson_vect
virtual Param * donne_para_poisson_vect(Param &para, int i) const =0
Function intended to be used by Map::poisson_vect and Map::poisson_vect_oohara .
Lorene::Tenseur::mp
const Map *const mp
Reference mapping.
Definition: tenseur.h:306
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::metric
const Metrique * metric
For tensor densities: the metric defining the conformal factor.
Definition: tenseur.h:325