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
24char 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// -----------------------
117namespace Lorene {
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 ...
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// -----------------------
201 Tenseur& scalaire) const {
202
203 Param bidon ;
205 resu.set_etat_qcq() ;
207 return resu ;
208}
209
210
211 //-----------------------------------//
212 // Vectorial Poisson equation //
213 // using Oohara scheme //
214 //-----------------------------------//
215
216// Version avec parametres
217// -----------------------
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// -----------------------
291
292 Param bidon ;
294 resu.set_etat_qcq() ;
296 return resu ;
297}
298
299
300 //---------------------------------------------//
301 // Vectorial Poisson equation TAU method//
302 //---------------------------------------------//
303
304// Version avec parametres
305// -----------------------
306void 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// -----------------------
388Tenseur 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// -----------------------
406void 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// -----------------------
476Tenseur 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}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Parameter storage.
Definition param.h:125
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
int get_type_indice(int i) const
Returns the type of the index number i .
Definition tenseur.h:726
const Metrique * metric
For tensor densities: the metric defining the conformal factor.
Definition tenseur.h:325
const Map *const mp
Reference mapping.
Definition tenseur.h:306
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
void poisson_vect_oohara(double lambda, Param &par, Tenseur &shift, Tenseur &scal) const
Solves the vectorial Poisson equation .
Tenseur(const Map &map, const Metrique *met=0x0, double weight=0)
Constructor for a scalar field.
Definition tenseur.C:209
friend Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tenseur.h:312
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition tenseur.C:1542
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition tenseur.C:645
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
int valence
Valence.
Definition tenseur.h:307
int etat
Logical state ETATZERO , ETATQCQ or ETATNONDEF.
Definition tenseur.h:321
double poids
For tensor densities: the weight.
Definition tenseur.h:323
int get_valence() const
Returns the valence.
Definition tenseur.h:710
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition tenseur.C:674
void poisson_vect(double lambda, Param &par, Tenseur &shift, Tenseur &vect, Tenseur &scal) const
Solves the vectorial Poisson equation : .
friend Tenseur skxk(const Tenseur &)
Contraction of the last index of (*this) with or , depending on the type of S .
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
Lorene prototypes.
Definition app_hor.h:64