LORENE
tenseur_pde_ylm.C
1/*
2 * Methods of the class tenseur for solving vectorial Poisson equations
3 * with a multipole falloff condition at the outer boundary
4 *
5 * (see file tenseur.h for documentation).
6 *
7 */
8
9/*
10 * Copyright (c) 2004 Joshua A. Faber
11 *
12 * This file is part of LORENE.
13 *
14 * LORENE is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU General Public License version 2
16 * as published by the Free Software Foundation.
17 *
18 * LORENE is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU General Public License for more details.
22 *
23 * You should have received a copy of the GNU General Public License
24 * along with LORENE; if not, write to the Free Software
25 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 *
27 */
28
29char tenseur_pde_ylm_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_pde_ylm.C,v 1.2 2014/10/13 08:53:42 j_novak Exp $" ;
30
31/*
32 * $Id: tenseur_pde_ylm.C,v 1.2 2014/10/13 08:53:42 j_novak Exp $
33 * $Log: tenseur_pde_ylm.C,v $
34 * Revision 1.2 2014/10/13 08:53:42 j_novak
35 * Lorene classes and functions now belong to the namespace Lorene.
36 *
37 * Revision 1.1 2004/12/29 16:32:33 k_taniguchi
38 * *** empty log message ***
39 *
40 *
41 * $Header: /cvsroot/Lorene/C++/Source/Tenseur/tenseur_pde_ylm.C,v 1.2 2014/10/13 08:53:42 j_novak Exp $
42 *
43 */
44
45// Lorene headers
46#include "map.h"
47#include "cmp.h"
48#include "param.h"
49#include "tenseur.h"
50
51 //-----------------------------------//
52 // Vectorial Poisson equation //
53 //-----------------------------------//
54
55// Version avec parametres
56// -----------------------
57namespace Lorene {
58void Tenseur::poisson_vect_ylm(double lambda, Param& para, Tenseur& shift,
59 Tenseur& vecteur, Tenseur& scalaire, int nylm,
60 double* intvec) const {
61 assert (lambda != -1) ;
62
63 // Verifications d'usage ...
64 assert (valence == 1) ;
65 assert (shift.get_valence() == 1) ;
66 assert (shift.get_type_indice(0) == type_indice(0)) ;
67 assert (vecteur.get_valence() == 1) ;
68 assert (vecteur.get_type_indice(0) == type_indice(0)) ;
69 assert (scalaire.get_valence() == 0) ;
70 assert (etat != ETATNONDEF) ;
71
72 Map_af mapping (*mp);
73
74 // Nothing to do if the source is zero
75 if (etat == ETATZERO) {
76
77 shift.set_etat_zero() ;
78
79 vecteur.set_etat_qcq() ;
80 for (int i=0; i<3; i++) {
81 vecteur.set(i) = 0 ;
82 }
83
84 scalaire.set_etat_qcq() ;
85 scalaire.set() = 0 ;
86
87 return ;
88 }
89
90 // On construit le tableau contenant le terme P_i ...
91 for (int i=0 ; i<3 ; i++) {
92 Param* par = mp->donne_para_poisson_vect(para, i) ;
93
94 double* intvec2=new double [nylm];
95 for (int j=0; j<nylm; j++) {
96 intvec2[j]=intvec[i*nylm+j];
97 }
98
99 (*this)(i).poisson_ylm(*par, vecteur.set(i),nylm,intvec2) ;
100
101 delete [] intvec2;
102
103 if (par != 0x0)
104 delete par ;
105 }
106 vecteur.set_triad( *triad ) ;
107
108 // Equation de Poisson scalaire :
109 Tenseur source_scal (-skxk(*this)) ;
110
111 Param* par = mp->donne_para_poisson_vect(para, 3) ;
112
113 double* intvec2=new double[nylm];
114 for (int j=0; j<nylm; j++) {
115 intvec2[j]=intvec[3*nylm+j];
116 }
117
118 source_scal().poisson_ylm(*par, scalaire.set(), nylm, intvec2) ;
119
120 delete [] intvec2;
121 if (par !=0x0)
122 delete par ;
123
124 // On construit le tableau contenant le terme d xsi / d x_i ...
125 Tenseur auxiliaire(scalaire) ;
126 Tenseur dxsi (auxiliaire.gradient()) ;
127
128 // On construit le tableau contenant le terme x_k d P_k / d x_i
129 Tenseur dp (skxk(vecteur.gradient())) ;
130
131 // Il ne reste plus qu'a tout ranger dans P :
132 // The final computation is done component by component because
133 // d_khi and x_d_w are covariant comp. whereas w_shift is
134 // contravariant
135
136 shift.set_etat_qcq() ;
137
138 for (int i=0 ; i<3 ; i++)
139 shift.set(i) = (lambda+2)/2/(lambda+1) * vecteur(i)
140 - (lambda/2/(lambda+1)) * (dxsi(i) + dp(i)) ;
141
142 shift.set_triad( *(vecteur.triad) ) ;
143
144}
145
146
147// Version sans parametres
148// -----------------------
149Tenseur Tenseur::poisson_vect_ylm(double lambda, Tenseur& vecteur,
150 Tenseur& scalaire, int nylm, double* intvec) const {
151
152 Param bidon ;
154 resu.set_etat_qcq() ;
155 poisson_vect_ylm(lambda, bidon, resu, vecteur, scalaire, nylm, intvec) ;
156 return resu ;
157}
158
159}
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
Tenseur(const Map &map, const Metrique *met=0x0, double weight=0)
Constructor for a scalar field.
Definition tenseur.C:209
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tenseur.h:312
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
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