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
32char 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// -----------------------
70namespace Lorene {
71void Tenseur::poisson_vect_regu(int k_div, int nzet, double unsgam1,
72 double lambda, Param& para, Tenseur& shift,
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
149 Tenseur dscal_div(*mp, 1, CON, mp->get_bvect_cart(), metric, poids) ;
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(),
164
165 delete par ;
166
167 // On construit le tableau contenant le terme d xsi / d x_i ...
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}
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
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
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tenseur.h:312
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_regu(int k_div, int nzet, double unsgam1, 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 .
Lorene prototypes.
Definition app_hor.h:64