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