LORENE
tenseur_pde.C
1 /*
2  * Copyright (c) 2000-2001 Eric Gourgoulhon
3  * Copyright (c) 2000-2001 Philippe Grandclement
4  *
5  * This file is part of LORENE.
6  *
7  * LORENE is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 2 of the License, or
10  * (at your option) any later version.
11  *
12  * LORENE is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with LORENE; if not, write to the Free Software
19  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20  *
21  */
22 
23 
24 
25 
26 /*
27  * $Id: tenseur_pde.C,v 1.8 2016/12/05 16:18:17 j_novak Exp $
28  * $Log: tenseur_pde.C,v $
29  * Revision 1.8 2016/12/05 16:18:17 j_novak
30  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
31  *
32  * Revision 1.7 2014/10/13 08:53:42 j_novak
33  * Lorene classes and functions now belong to the namespace Lorene.
34  *
35  * Revision 1.6 2006/06/01 12:47:54 p_grandclement
36  * update of the Bin_ns_bh project
37  *
38  * Revision 1.5 2005/08/30 08:35:13 p_grandclement
39  * Addition of the Tau version of the vectorial Poisson equation for the Tensors
40  *
41  * Revision 1.4 2003/10/03 15:58:51 j_novak
42  * Cleaning of some headers
43  *
44  * Revision 1.3 2002/08/07 16:14:11 j_novak
45  * class Tenseur can now also handle tensor densities, this should be transparent to older codes
46  *
47  * Revision 1.2 2002/07/09 16:46:23 p_grandclement
48  * The Param in the case of an affine mapping is now 0x0 and not deleted
49  * (I wonder why it was working before)
50  *
51  * Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
52  * LORENE
53  *
54  * Revision 2.13 2000/10/04 14:58:32 eric
55  * Ajout de shift.set_etat_qcq() avant l'affection de shift.
56  *
57  * Revision 2.12 2000/09/27 15:07:22 eric
58  * Traitement source nulle dans poisson_vect.
59  *
60  * Revision 2.11 2000/05/22 15:48:18 phil
61  * modification de Oohara pour passer avec dzpuis == 2
62  * pardon 3
63  *
64  * Revision 2.10 2000/05/22 15:00:46 phil
65  * Modification de la methode de Shibata :
66  * on doit passer une source en r^4 et l equation scalaire est alors resolue en utilisant l'algotihme avec dzpuis == 3
67  *
68  * Revision 2.9 2000/03/10 15:51:42 eric
69  * Traitement dzpuis de source_scal.
70  *
71  * Revision 2.8 2000/03/08 10:21:05 eric
72  * Appel de delete sur le Param* retourne par Map_et::donne_para_poisson_vect[
73  * lorsqu'il n'est plus utilise (correction Memory leak).
74  *
75  * Revision 2.7 2000/03/07 16:53:42 eric
76  * *** empty log message ***
77  *
78  * Revision 2.6 2000/03/07 15:43:32 phil
79  * gestion des cas dzpuis ==4
80  *
81  * Revision 2.5 2000/02/21 12:55:09 eric
82  * Traitement des triades.
83  *
84  * Revision 2.4 2000/02/16 17:13:05 eric
85  * Correction
86  * mp->donne_para_poisson_vect(para, 4)
87  * devient
88  * mp->donne_para_poisson_vect(para, 3)
89  * dans Tenseur::poisson_vect
90  * /
91  *
92  * Revision 2.3 2000/02/15 10:26:49 phil
93  * le calcul de sol n'appelle plus Map::poisson_vect mais est fait dans
94  * Tenseur::poisson_vect (respectivement poisson_vect_oohara)
95  *
96  * Revision 2.2 2000/02/09 19:32:58 eric
97  * La triade de decomposition est desormais passee en argument des constructeurs.
98  *
99  * Revision 2.1 2000/02/09 10:01:32 phil
100  * ajout version oohara
101  *
102  * Revision 2.0 2000/01/21 12:58:57 phil
103  * *** empty log message ***
104  *
105  *
106  * $Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_pde.C,v 1.8 2016/12/05 16:18:17 j_novak Exp $
107  *
108  */
109 
110 // Header Lorene:
111 #include "param.h"
112 #include "tenseur.h"
113 
114  //-----------------------------------//
115  // Vectorial Poisson equation //
116  //-----------------------------------//
117 
118 // Version avec parametres
119 // -----------------------
120 namespace Lorene {
121 void Tenseur::poisson_vect(double lambda, Param& para, Tenseur& shift
122  , Tenseur& vecteur, Tenseur& scalaire) const {
123  assert (lambda != -1) ;
124 
125  // Verifications d'usage ...
126  assert (valence == 1) ;
127  assert (shift.get_valence() == 1) ;
128  assert (shift.get_type_indice(0) == type_indice(0)) ;
129  assert (vecteur.get_valence() == 1) ;
130  assert (vecteur.get_type_indice(0) == type_indice(0)) ;
131  assert (scalaire.get_valence() == 0) ;
132  assert (etat != ETATNONDEF) ;
133 
134  // Nothing to do if the source is zero
135  if (etat == ETATZERO) {
136 
137  shift.set_etat_zero() ;
138 
139  vecteur.set_etat_qcq() ;
140  for (int i=0; i<3; i++) {
141  vecteur.set(i) = 0 ;
142  }
143 
144  scalaire.set_etat_qcq() ;
145  scalaire.set() = 0 ;
146 
147  return ;
148  }
149 
150  for (int i=0 ; i<3 ; i++)
151  assert ((*this)(i).check_dzpuis(4)) ;
152 
153  // On construit le tableau contenant le terme P_i ...
154  for (int i=0 ; i<3 ; i++) {
155  Param* par = mp->donne_para_poisson_vect(para, i) ;
156 
157  (*this)(i).poisson(*par, vecteur.set(i)) ;
158 
159  if (par != 0x0)
160  delete par ;
161  }
162  vecteur.set_triad( *triad ) ;
163 
164  // Equation de Poisson scalaire :
165  Tenseur source_scal (-skxk(*this)) ;
166 
167  assert (source_scal().check_dzpuis(3)) ;
168 
169  Param* par = mp->donne_para_poisson_vect(para, 3) ;
170 
171  source_scal().poisson(*par, scalaire.set()) ;
172 
173  if (par !=0x0)
174  delete par ;
175 
176  // On construit le tableau contenant le terme d xsi / d x_i ...
177  Tenseur auxiliaire(scalaire) ;
178  Tenseur dxsi (auxiliaire.gradient()) ;
179  dxsi.dec2_dzpuis() ;
180 
181  // On construit le tableau contenant le terme x_k d P_k / d x_i
182  Tenseur dp (skxk(vecteur.gradient())) ;
183  dp.dec_dzpuis() ;
184 
185  // Il ne reste plus qu'a tout ranger dans P :
186  // The final computation is done component by component because
187  // d_khi and x_d_w are covariant comp. whereas w_shift is
188  // contravariant
189 
190  shift.set_etat_qcq() ;
191 
192  for (int i=0 ; i<3 ; i++)
193  shift.set(i) = (lambda+2)/2/(lambda+1) * vecteur(i)
194  - (lambda/2/(lambda+1)) * (dxsi(i) + dp(i)) ;
195 
196  shift.set_triad( *(vecteur.triad) ) ;
197 
198 }
199 
200 
201 // Version sans parametres
202 // -----------------------
203 Tenseur Tenseur::poisson_vect(double lambda, Tenseur& vecteur,
204  Tenseur& scalaire) const {
205 
206  Param bidon ;
208  resu.set_etat_qcq() ;
209  poisson_vect(lambda, bidon, resu, vecteur, scalaire) ;
210  return resu ;
211 }
212 
213 
214  //-----------------------------------//
215  // Vectorial Poisson equation //
216  // using Oohara scheme //
217  //-----------------------------------//
218 
219 // Version avec parametres
220 // -----------------------
221 void Tenseur::poisson_vect_oohara(double lambda, Param& para, Tenseur& shift,
222  Tenseur& chi) const {
223 
224  // Ne marche pas pour lambda =-1
225  assert (lambda != -1) ;
226 
227  // Verifications d'usage ...
228  assert (valence == 1) ;
229  assert (shift.get_valence() == 1) ;
230  assert (shift.get_type_indice(0) == type_indice(0)) ;
231  assert (chi.get_valence() == 0) ;
232  assert (etat != ETATNONDEF) ;
233 
234  // Nothing to do if the source is zero
235  if (etat == ETATZERO) {
236  shift.set_etat_zero() ;
237  chi.set_etat_qcq() ;
238  chi.set() = 0 ;
239  return ;
240  }
241 
242  for (int i=0 ; i<3 ; i++)
243  assert ((*this)(i).check_dzpuis(3) ||
244  (*this)(i).check_dzpuis(4)) ;
245 
246 
247  Tenseur copie(*this) ;
248  copie.dec2_dzpuis() ;
249  if ((*this)(0).check_dzpuis(4))
250  copie.dec2_dzpuis() ;
251  else
252  copie.dec_dzpuis() ;
253 
254  Tenseur source_scal(contract(copie.gradient(), 0, 1)/(1.+lambda)) ;
255  source_scal.inc2_dzpuis() ;
256 
257  Param* par = mp->donne_para_poisson_vect(para, 3) ;
258 
259  source_scal().poisson(*par, chi.set());
260  if (par !=0x0)
261  delete par ;
262 
263  Tenseur source_vect(*this) ;
264  if ((*this)(0).check_dzpuis(4))
265  source_vect.dec_dzpuis() ;
266  Tenseur chi_grad (chi.gradient()) ;
267  chi_grad.inc_dzpuis() ;
268 
269  for (int i=0 ; i<3 ; i++)
270  source_vect.set(i) -= lambda*chi_grad(i) ;
271  assert( *(source_vect.triad) == *((chi.gradient()).get_triad()) ) ;
272 
273  if (shift.get_etat() == ETATZERO) {
274  shift.set_etat_qcq() ;
275  for (int i=0 ; i<3 ; i++)
276  shift.set(i) = 0 ;
277  }
278 
279  for (int i=0 ; i<3 ; i++) {
280  par = mp->donne_para_poisson_vect(para, i) ;
281  source_vect(i).poisson(*par, shift.set(i)) ;
282 
283  if (par !=0x0)
284  delete par ;
285  }
286  shift.set_triad( *(source_vect.triad) ) ;
287 
288 }
289 
290 
291 // Version sans parametres
292 // -----------------------
293 Tenseur Tenseur::poisson_vect_oohara(double lambda, Tenseur& scalaire) const {
294 
295  Param bidon ;
297  resu.set_etat_qcq() ;
298  poisson_vect_oohara(lambda, bidon, resu, scalaire) ;
299  return resu ;
300 }
301 
302 
303  //---------------------------------------------//
304  // Vectorial Poisson equation TAU method//
305  //---------------------------------------------//
306 
307 // Version avec parametres
308 // -----------------------
309 void Tenseur::poisson_vect_tau(double lambda, Param& para, Tenseur& shift
310  , Tenseur& vecteur, Tenseur& scalaire) const {
311  assert (lambda != -1) ;
312 
313  // Verifications d'usage ...
314  assert (valence == 1) ;
315  assert (shift.get_valence() == 1) ;
316  assert (shift.get_type_indice(0) == type_indice(0)) ;
317  assert (vecteur.get_valence() == 1) ;
318  assert (vecteur.get_type_indice(0) == type_indice(0)) ;
319  assert (scalaire.get_valence() == 0) ;
320  assert (etat != ETATNONDEF) ;
321 
322  // Nothing to do if the source is zero
323  if (etat == ETATZERO) {
324 
325  shift.set_etat_zero() ;
326 
327  vecteur.set_etat_qcq() ;
328  for (int i=0; i<3; i++) {
329  vecteur.set(i) = 0 ;
330  }
331 
332  scalaire.set_etat_qcq() ;
333  scalaire.set() = 0 ;
334 
335  return ;
336  }
337 
338  for (int i=0 ; i<3 ; i++)
339  assert ((*this)(i).check_dzpuis(4)) ;
340 
341  // On construit le tableau contenant le terme P_i ...
342  for (int i=0 ; i<3 ; i++) {
343  Param* par = mp->donne_para_poisson_vect(para, i) ;
344 
345  (*this)(i).poisson_tau(*par, vecteur.set(i)) ;
346 
347  if (par != 0x0)
348  delete par ;
349  }
350  vecteur.set_triad( *triad ) ;
351 
352  // Equation de Poisson scalaire :
353  Tenseur source_scal (-skxk(*this)) ;
354 
355  assert (source_scal().check_dzpuis(3)) ;
356 
357  Param* par = mp->donne_para_poisson_vect(para, 3) ;
358 
359  source_scal().poisson_tau(*par, scalaire.set()) ;
360 
361  if (par !=0x0)
362  delete par ;
363 
364  // On construit le tableau contenant le terme d xsi / d x_i ...
365  Tenseur auxiliaire(scalaire) ;
366  Tenseur dxsi (auxiliaire.gradient()) ;
367  dxsi.dec2_dzpuis() ;
368 
369  // On construit le tableau contenant le terme x_k d P_k / d x_i
370  Tenseur dp (skxk(vecteur.gradient())) ;
371  dp.dec_dzpuis() ;
372 
373  // Il ne reste plus qu'a tout ranger dans P :
374  // The final computation is done component by component because
375  // d_khi and x_d_w are covariant comp. whereas w_shift is
376  // contravariant
377 
378  shift.set_etat_qcq() ;
379 
380  for (int i=0 ; i<3 ; i++)
381  shift.set(i) = (lambda+2)/2/(lambda+1) * vecteur(i)
382  - (lambda/2/(lambda+1)) * (dxsi(i) + dp(i)) ;
383 
384  shift.set_triad( *(vecteur.triad) ) ;
385 
386 }
387 
388 
389 // Version sans parametres
390 // -----------------------
391 Tenseur Tenseur::poisson_vect_tau(double lambda, Tenseur& vecteur,
392  Tenseur& scalaire) const {
393 
394  Param bidon ;
396  resu.set_etat_qcq() ;
397  poisson_vect_tau(lambda, bidon, resu, vecteur, scalaire) ;
398  return resu ;
399 }
400 
401 
402  //-----------------------------------//
403  // Vectorial Poisson equation //
404  // using Oohara scheme //
405  //-----------------------------------//
406 
407 // Version avec parametres
408 // -----------------------
409 void Tenseur::poisson_vect_oohara_tau(double lambda, Param& para, Tenseur& shift,
410  Tenseur& chi) const {
411 
412  // Ne marche pas pour lambda =-1
413  assert (lambda != -1) ;
414 
415  // Verifications d'usage ...
416  assert (valence == 1) ;
417  assert (shift.get_valence() == 1) ;
418  assert (shift.get_type_indice(0) == type_indice(0)) ;
419  assert (chi.get_valence() == 0) ;
420  assert (etat != ETATNONDEF) ;
421 
422  // Nothing to do if the source is zero
423  if (etat == ETATZERO) {
424  shift.set_etat_zero() ;
425  chi.set_etat_qcq() ;
426  chi.set() = 0 ;
427  return ;
428  }
429 
430  for (int i=0 ; i<3 ; i++)
431  assert ((*this)(i).check_dzpuis(3) ||
432  (*this)(i).check_dzpuis(4)) ;
433 
434 
435  Tenseur copie(*this) ;
436  copie.dec2_dzpuis() ;
437  if ((*this)(0).check_dzpuis(4))
438  copie.dec2_dzpuis() ;
439  else
440  copie.dec_dzpuis() ;
441 
442  Tenseur source_scal(contract(copie.gradient(), 0, 1)/(1.+lambda)) ;
443  source_scal.inc2_dzpuis() ;
444 
445  Param* par = mp->donne_para_poisson_vect(para, 3) ;
446 
447  source_scal().poisson_tau(*par, chi.set());
448 
449  if (par !=0x0)
450  delete par ;
451 
452  Tenseur source_vect(*this) ;
453  if ((*this)(0).check_dzpuis(4))
454  source_vect.dec_dzpuis() ;
455 
456  Tenseur chi_grad (chi.gradient()) ;
457  chi_grad.inc_dzpuis() ;
458 
459  for (int i=0 ; i<3 ; i++)
460  source_vect.set(i) -= lambda*chi_grad(i) ;
461 
462  assert( *(source_vect.triad) == *((chi.gradient()).get_triad()) ) ;
463 
464  for (int i=0 ; i<3 ; i++) {
465  par = mp->donne_para_poisson_vect(para, i) ;
466 
467  source_vect(i).poisson_tau(*par, shift.set(i)) ;
468 
469  if (par !=0x0)
470  delete par ;
471  }
472  shift.set_triad( *(source_vect.triad) ) ;
473 
474 }
475 
476 
477 // Version sans parametres
478 // -----------------------
479 Tenseur Tenseur::poisson_vect_oohara_tau(double lambda, Tenseur& scalaire) const {
480 
481  Param bidon ;
483  resu.set_etat_qcq() ;
484  poisson_vect_oohara_tau(lambda, bidon, resu, scalaire) ;
485  return resu ;
486 }
487 
488 }
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
void dec2_dzpuis()
dzpuis -= 2 ;
Definition: tenseur.C:1146
void poisson_vect(double lambda, Param &par, Tenseur &shift, Tenseur &vect, Tenseur &scal) const
Solves the vectorial Poisson equation : .
Definition: tenseur_pde.C:121
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition: tenseur.C:690
void poisson_vect_oohara(double lambda, Param &par, Tenseur &shift, Tenseur &scal) const
Solves the vectorial Poisson equation .
Definition: tenseur_pde.C:221
Lorene prototypes.
Definition: app_hor.h:67
int etat
Logical state ETATZERO , ETATQCQ or ETATNONDEF.
Definition: tenseur.h:324
void inc_dzpuis()
dzpuis += 1 ;
Definition: tenseur.C:1133
void inc2_dzpuis()
dzpuis += 2 ;
Definition: tenseur.C:1159
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
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tenseur.h:707
double poids
For tensor densities: the weight.
Definition: tenseur.h:326
Parameter storage.
Definition: param.h:125
int get_etat() const
Returns the logical state.
Definition: tenseur.h:710
friend Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
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 .
Tenseur(const Map &map, const Metrique *met=0x0, double weight=0)
Constructor for a scalar field.
Definition: tenseur.C:225
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
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