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 
33 
34 /*
35  * $Id: tenseur_pde_regu.C,v 1.5 2016/12/05 16:18:17 j_novak Exp $
36  * $Log: tenseur_pde_regu.C,v $
37  * Revision 1.5 2016/12/05 16:18:17 j_novak
38  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
39  *
40  * Revision 1.4 2014/10/13 08:53:42 j_novak
41  * Lorene classes and functions now belong to the namespace Lorene.
42  *
43  * Revision 1.3 2003/10/03 15:58:51 j_novak
44  * Cleaning of some headers
45  *
46  * Revision 1.2 2002/08/07 16:14:11 j_novak
47  * class Tenseur can now also handle tensor densities, this should be transparent to older codes
48  *
49  * Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
50  * LORENE
51  *
52  * Revision 2.1 2001/01/15 11:01:34 phil
53  * vire version sans parametres
54  *
55  * Revision 2.0 2000/10/06 15:34:03 keisuke
56  * Initial revision.
57  *
58  *
59  * $Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_pde_regu.C,v 1.5 2016/12/05 16:18:17 j_novak Exp $
60  *
61  */
62 
63 // Header Lorene:
64 #include "param.h"
65 #include "tenseur.h"
66 
67  //-----------------------------------//
68  // Vectorial Poisson equation //
69  //-----------------------------------//
70 
71 // Version avec parametres
72 // -----------------------
73 namespace Lorene {
74 void Tenseur::poisson_vect_regu(int k_div, int nzet, double unsgam1,
75  double lambda, Param& para, Tenseur& shift,
76  Tenseur& vecteur, Tenseur& scalaire) const {
77  assert (lambda != -1) ;
78 
79  // Verifications d'usage ...
80  assert (valence == 1) ;
81  assert (shift.get_valence() == 1) ;
82  assert (shift.get_type_indice(0) == type_indice(0)) ;
83  assert (vecteur.get_valence() == 1) ;
84  assert (vecteur.get_type_indice(0) == type_indice(0)) ;
85  assert (scalaire.get_valence() == 0) ;
86  assert (etat != ETATNONDEF) ;
87 
88  // Nothing to do if the source is zero
89  if (etat == ETATZERO) {
90 
91  shift.set_etat_zero() ;
92 
93  vecteur.set_etat_qcq() ;
94  for (int i=0; i<3; i++) {
95  vecteur.set(i) = 0 ;
96  }
97 
98  scalaire.set_etat_qcq() ;
99  scalaire.set() = 0 ;
100 
101  return ;
102  }
103 
104  for (int i=0 ; i<3 ; i++)
105  assert ((*this)(i).check_dzpuis(4)) ;
106 
107  Tenseur vecteur_regu(*mp, 1, CON, mp->get_bvect_cart(), metric, poids) ;
108  Tenseur vecteur_div(*mp, 1, CON, mp->get_bvect_cart(), metric, poids) ;
109  Tenseur dvect_div(*mp, 1, CON, mp->get_bvect_cart(), metric, poids) ;
110  Tenseur souvect_regu(*mp, 1, CON, mp->get_bvect_cart(), metric, poids) ;
111  Tenseur souvect_div(*mp, 1, CON, mp->get_bvect_cart(), metric, poids) ;
112 
113  vecteur_regu.set_etat_qcq() ;
114  vecteur_div.set_etat_qcq() ;
115  dvect_div.set_etat_qcq() ;
116  souvect_regu.set_etat_qcq() ;
117  souvect_div.set_etat_qcq() ;
118 
119  // On construit le tableau contenant le terme P_i ...
120 
121  // Apply only to x and y components because poisson_regular does not
122  // work for z component due to the symmetry.
123  for (int i=0 ; i<2 ; i++) {
124  Param* par = mp->donne_para_poisson_vect(para, i) ;
125 
126  (*this)(i).poisson_regular(k_div, nzet, unsgam1, *par,
127  vecteur.set(i),
128  vecteur_regu.set(i), vecteur_div.set(i),
129  dvect_div,
130  souvect_regu.set(i), souvect_div.set(i)) ;
131 
132  delete par ;
133  }
134 
135  Param* par = mp->donne_para_poisson_vect(para, 2) ;
136 
137  (*this)(2).poisson(*par, vecteur.set(2)) ;
138 
139  delete par ;
140 
141  vecteur.set_triad( *triad ) ;
142 
143  // Equation de Poisson scalaire :
144  Tenseur source_scal (-skxk(*this)) ;
145 
146  assert (source_scal().check_dzpuis(3)) ;
147 
148  par = mp->donne_para_poisson_vect(para, 3) ;
149 
150  Tenseur scalaire_regu(*mp, metric, poids) ;
151  Tenseur scalaire_div(*mp, metric, poids) ;
152  Tenseur dscal_div(*mp, 1, CON, mp->get_bvect_cart(), metric, poids) ;
153  Cmp souscal_regu(mp) ;
154  Cmp souscal_div(mp) ;
155 
156  scalaire_regu.set_etat_qcq() ;
157  scalaire_div.set_etat_qcq() ;
158  dscal_div.set_etat_qcq() ;
159 
160  souscal_regu.std_base_scal() ;
161  souscal_div.std_base_scal() ;
162 
163  source_scal().poisson_regular(k_div, nzet, unsgam1, *par,
164  scalaire.set(),
165  scalaire_regu.set(), scalaire_div.set(),
166  dscal_div, souscal_regu, souscal_div) ;
167 
168  delete par ;
169 
170  // On construit le tableau contenant le terme d xsi / d x_i ...
171  Tenseur auxiliaire(scalaire) ;
172  Tenseur dxsi (auxiliaire.gradient()) ;
173  dxsi.dec2_dzpuis() ;
174 
175  // On construit le tableau contenant le terme x_k d P_k / d x_i
176  Tenseur dp (skxk(vecteur.gradient())) ;
177  dp.dec_dzpuis() ;
178 
179  // Il ne reste plus qu'a tout ranger dans P :
180  // The final computation is done component by component because
181  // d_khi and x_d_w are covariant comp. whereas w_shift is
182  // contravariant
183 
184  shift.set_etat_qcq() ;
185 
186  for (int i=0 ; i<3 ; i++)
187  shift.set(i) = (lambda+2)/2/(lambda+1) * vecteur(i)
188  - (lambda/2/(lambda+1)) * (dxsi(i) + dp(i)) ;
189 
190  shift.set_triad( *(vecteur.triad) ) ;
191 
192 }
193 
194 
195 }
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:321
const Map *const mp
Reference mapping.
Definition: tenseur.h:309
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tenseur.h:315
int get_type_indice(int i) const
Returns the type of the index number i .
Definition: tenseur.h:729
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
void dec2_dzpuis()
dzpuis -= 2 ;
Definition: tenseur.C:1146
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition: tenseur.C:690
Lorene prototypes.
Definition: app_hor.h:67
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
int etat
Logical state ETATZERO , ETATQCQ or ETATNONDEF.
Definition: tenseur.h:324
int get_valence() const
Returns the valence.
Definition: tenseur.h:713
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:840
int valence
Valence.
Definition: tenseur.h:310
double poids
For tensor densities: the weight.
Definition: tenseur.h:326
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:647
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:724
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:809
void dec_dzpuis()
dzpuis -= 1 ;
Definition: tenseur.C:1120
const Metrique * metric
For tensor densities: the metric defining the conformal factor.
Definition: tenseur.h:328
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:652
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition: tenseur.C:661
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:304
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition: tenseur.C:1558