LORENE
tenseur_pde_regu.C
1 /*
2  * Method of class Tenseur to solve a vector Poisson equation
3  * by regularizing its source.
4  *
5  * (see file tenseur.h for documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 2000-2001 Keisuke Taniguchi
11  * Copyright (c) 2001 Philippe Grandclement
12  *
13  * This file is part of LORENE.
14  *
15  * LORENE is free software; you can redistribute it and/or modify
16  * it under the terms of the GNU General Public License as published by
17  * the Free Software Foundation; either version 2 of the License, or
18  * (at your option) any later version.
19  *
20  * LORENE is distributed in the hope that it will be useful,
21  * but WITHOUT ANY WARRANTY; without even the implied warranty of
22  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23  * GNU General Public License for more details.
24  *
25  * You should have received a copy of the GNU General Public License
26  * along with LORENE; if not, write to the Free Software
27  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28  *
29  */
30 
31 
32 char tenseur_pde_regu_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_pde_regu.C,v 1.4 2014/10/13 08:53:42 j_novak Exp $" ;
33 
34 /*
35  * $Id: tenseur_pde_regu.C,v 1.4 2014/10/13 08:53:42 j_novak Exp $
36  * $Log: tenseur_pde_regu.C,v $
37  * Revision 1.4 2014/10/13 08:53:42 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.3 2003/10/03 15:58:51 j_novak
41  * Cleaning of some headers
42  *
43  * Revision 1.2 2002/08/07 16:14:11 j_novak
44  * class Tenseur can now also handle tensor densities, this should be transparent to older codes
45  *
46  * Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
47  * LORENE
48  *
49  * Revision 2.1 2001/01/15 11:01:34 phil
50  * vire version sans parametres
51  *
52  * Revision 2.0 2000/10/06 15:34:03 keisuke
53  * Initial revision.
54  *
55  *
56  * $Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_pde_regu.C,v 1.4 2014/10/13 08:53:42 j_novak Exp $
57  *
58  */
59 
60 // Header Lorene:
61 #include "param.h"
62 #include "tenseur.h"
63 
64  //-----------------------------------//
65  // Vectorial Poisson equation //
66  //-----------------------------------//
67 
68 // Version avec parametres
69 // -----------------------
70 namespace Lorene {
71 void Tenseur::poisson_vect_regu(int k_div, int nzet, double unsgam1,
72  double lambda, Param& para, Tenseur& shift,
73  Tenseur& vecteur, Tenseur& scalaire) const {
74  assert (lambda != -1) ;
75 
76  // Verifications d'usage ...
77  assert (valence == 1) ;
78  assert (shift.get_valence() == 1) ;
79  assert (shift.get_type_indice(0) == type_indice(0)) ;
80  assert (vecteur.get_valence() == 1) ;
81  assert (vecteur.get_type_indice(0) == type_indice(0)) ;
82  assert (scalaire.get_valence() == 0) ;
83  assert (etat != ETATNONDEF) ;
84 
85  // Nothing to do if the source is zero
86  if (etat == ETATZERO) {
87 
88  shift.set_etat_zero() ;
89 
90  vecteur.set_etat_qcq() ;
91  for (int i=0; i<3; i++) {
92  vecteur.set(i) = 0 ;
93  }
94 
95  scalaire.set_etat_qcq() ;
96  scalaire.set() = 0 ;
97 
98  return ;
99  }
100 
101  for (int i=0 ; i<3 ; i++)
102  assert ((*this)(i).check_dzpuis(4)) ;
103 
104  Tenseur vecteur_regu(*mp, 1, CON, mp->get_bvect_cart(), metric, poids) ;
105  Tenseur vecteur_div(*mp, 1, CON, mp->get_bvect_cart(), metric, poids) ;
106  Tenseur dvect_div(*mp, 1, CON, mp->get_bvect_cart(), metric, poids) ;
107  Tenseur souvect_regu(*mp, 1, CON, mp->get_bvect_cart(), metric, poids) ;
108  Tenseur souvect_div(*mp, 1, CON, mp->get_bvect_cart(), metric, poids) ;
109 
110  vecteur_regu.set_etat_qcq() ;
111  vecteur_div.set_etat_qcq() ;
112  dvect_div.set_etat_qcq() ;
113  souvect_regu.set_etat_qcq() ;
114  souvect_div.set_etat_qcq() ;
115 
116  // On construit le tableau contenant le terme P_i ...
117 
118  // Apply only to x and y components because poisson_regular does not
119  // work for z component due to the symmetry.
120  for (int i=0 ; i<2 ; i++) {
121  Param* par = mp->donne_para_poisson_vect(para, i) ;
122 
123  (*this)(i).poisson_regular(k_div, nzet, unsgam1, *par,
124  vecteur.set(i),
125  vecteur_regu.set(i), vecteur_div.set(i),
126  dvect_div,
127  souvect_regu.set(i), souvect_div.set(i)) ;
128 
129  delete par ;
130  }
131 
132  Param* par = mp->donne_para_poisson_vect(para, 2) ;
133 
134  (*this)(2).poisson(*par, vecteur.set(2)) ;
135 
136  delete par ;
137 
138  vecteur.set_triad( *triad ) ;
139 
140  // Equation de Poisson scalaire :
141  Tenseur source_scal (-skxk(*this)) ;
142 
143  assert (source_scal().check_dzpuis(3)) ;
144 
145  par = mp->donne_para_poisson_vect(para, 3) ;
146 
147  Tenseur scalaire_regu(*mp, metric, poids) ;
148  Tenseur scalaire_div(*mp, metric, poids) ;
149  Tenseur dscal_div(*mp, 1, CON, mp->get_bvect_cart(), metric, poids) ;
150  Cmp souscal_regu(mp) ;
151  Cmp souscal_div(mp) ;
152 
153  scalaire_regu.set_etat_qcq() ;
154  scalaire_div.set_etat_qcq() ;
155  dscal_div.set_etat_qcq() ;
156 
157  souscal_regu.std_base_scal() ;
158  souscal_div.std_base_scal() ;
159 
160  source_scal().poisson_regular(k_div, nzet, unsgam1, *par,
161  scalaire.set(),
162  scalaire_regu.set(), scalaire_div.set(),
163  dscal_div, souscal_regu, souscal_div) ;
164 
165  delete par ;
166 
167  // On construit le tableau contenant le terme d xsi / d x_i ...
168  Tenseur auxiliaire(scalaire) ;
169  Tenseur dxsi (auxiliaire.gradient()) ;
170  dxsi.dec2_dzpuis() ;
171 
172  // On construit le tableau contenant le terme x_k d P_k / d x_i
173  Tenseur dp (skxk(vecteur.gradient())) ;
174  dp.dec_dzpuis() ;
175 
176  // Il ne reste plus qu'a tout ranger dans P :
177  // The final computation is done component by component because
178  // d_khi and x_d_w are covariant comp. whereas w_shift is
179  // contravariant
180 
181  shift.set_etat_qcq() ;
182 
183  for (int i=0 ; i<3 ; i++)
184  shift.set(i) = (lambda+2)/2/(lambda+1) * vecteur(i)
185  - (lambda/2/(lambda+1)) * (dxsi(i) + dp(i)) ;
186 
187  shift.set_triad( *(vecteur.triad) ) ;
188 
189 }
190 
191 
192 }
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:315
const Map *const mp
Reference mapping.
Definition: tenseur.h:303
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tenseur.h:309
int get_type_indice(int i) const
Returns the type of the index number i .
Definition: tenseur.h:723
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:443
void dec2_dzpuis()
dzpuis -= 2 ;
Definition: tenseur.C:1130
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition: tenseur.C:674
Lorene prototypes.
Definition: app_hor.h:64
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
int etat
Logical state ETATZERO , ETATQCQ or ETATNONDEF.
Definition: tenseur.h:318
int get_valence() const
Returns the valence.
Definition: tenseur.h:707
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
int valence
Valence.
Definition: tenseur.h:304
double poids
For tensor densities: the weight.
Definition: tenseur.h:320
Parameter storage.
Definition: param.h:125
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 .
friend Tenseur skxk(const Tenseur &)
Contraction of the last index of (*this) with or , depending on the type of S .
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:644
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:721
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:791
void dec_dzpuis()
dzpuis -= 1 ;
Definition: tenseur.C:1104
const Metrique * metric
For tensor densities: the metric defining the conformal factor.
Definition: tenseur.h:322
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition: tenseur.C:645
void poisson_vect_regu(int k_div, int nzet, double unsgam1, double lambda, Param &par, Tenseur &shift, Tenseur &vect, Tenseur &scal) const
Solves the vectorial Poisson equation : .
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:298
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition: tenseur.C:1542