LORENE
sym_tensor_trans_pde.C
1/*
2 * Functions to solve various PDEs for a divergence-free symmetric tensor.
3 *
4 * (see file sym_tensor.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2005-2006 Jerome Novak
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28char sym_tensor_trans_pde_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_trans_pde.C,v 1.16 2014/10/13 08:53:43 j_novak Exp $" ;
29
30/*
31 * $Id: sym_tensor_trans_pde.C,v 1.16 2014/10/13 08:53:43 j_novak Exp $
32 * $Log: sym_tensor_trans_pde.C,v $
33 * Revision 1.16 2014/10/13 08:53:43 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.15 2014/10/06 15:13:19 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.14 2010/10/11 10:38:34 j_novak
40 * *** empty log message ***
41 *
42 * Revision 1.13 2010/10/11 10:23:03 j_novak
43 * Removed methods Sym_tensor_trans::solve_hrr() and Sym_tensor_trans::set_WX_det_one(), as they are no longer relevant.
44 *
45 * Revision 1.12 2006/09/05 15:38:45 j_novak
46 * The fuctions sol_Dirac... are in a seperate file, with new parameters to
47 * control the boundary conditions.
48 *
49 * Revision 1.11 2006/08/31 12:13:22 j_novak
50 * Added an argument of type Param to Sym_tensor_trans::sol_Dirac_A().
51 *
52 * Revision 1.10 2006/06/28 07:48:26 j_novak
53 * Better treatment of some null cases.
54 *
55 * Revision 1.9 2006/06/21 15:42:47 j_novak
56 * Minor changes.
57 *
58 * Revision 1.8 2006/06/20 12:07:15 j_novak
59 * Improved execution speed for sol_Dirac_tildeB...
60 *
61 * Revision 1.7 2006/06/14 10:04:21 j_novak
62 * New methods sol_Dirac_l01, set_AtB_det_one and set_AtB_trace_zero.
63 *
64 * Revision 1.6 2006/06/13 13:30:12 j_novak
65 * New members sol_Dirac_A and sol_Dirac_tildeB (see documentation).
66 *
67 * Revision 1.5 2006/06/12 13:37:23 j_novak
68 * Added bounds in l (multipolar momentum) for Sym_tensor_trans::solve_hrr.
69 *
70 * Revision 1.4 2005/11/28 14:45:17 j_novak
71 * Improved solution of the Poisson tensor equation in the case of a transverse
72 * tensor.
73 *
74 * Revision 1.3 2005/11/24 14:07:54 j_novak
75 * Use of Matrice::annule_hard()
76 *
77 * Revision 1.2 2005/11/24 09:24:25 j_novak
78 * Corrected some missing references.
79 *
80 * Revision 1.1 2005/09/16 13:58:11 j_novak
81 * New Poisson solver for a Sym_tensor_trans.
82 *
83 *
84 * $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_trans_pde.C,v 1.16 2014/10/13 08:53:43 j_novak Exp $
85 *
86 */
87
88// C headers
89#include <cassert>
90#include <cmath>
91
92// Lorene headers
93#include "tensor.h"
94#include "diff.h"
95#include "proto.h"
96#include "param.h"
97
98namespace Lorene {
100
101 // All this has a meaning only for spherical components...
102 assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
103 //## ... and affine mapping, for the moment!
104 const Map_af* mpaff = dynamic_cast<const Map_af*>(mp) ;
105 assert( mpaff!= 0x0) ;
107
108 const Mg3d& gri = *mp->get_mg() ;
109 int np = gri.get_np(0) ;
110 int nt = gri.get_nt(0) ;
111 assert (nt > 4) ;
112 if (np == 1) {
113 int nz = gri.get_nzone() ;
114 double* bornes = new double[nz+1] ;
115 const double* alp = mpaff->get_alpha() ;
116 const double* bet = mpaff->get_beta() ;
117 for (int lz=0; lz<nz; lz++) {
118 assert (gri.get_np(lz) == np) ;
119 assert (gri.get_nt(lz) == nt) ;
120 switch (gri.get_type_r(lz)) {
121 case RARE: {
122 bornes[lz] = bet[lz] ;
123 break ;
124 }
125 case FIN: {
126 bornes[lz] = bet[lz] - alp[lz] ;
127 break ;
128 }
129 case UNSURR: {
130 bornes[lz] = double(1) / ( bet[lz] - alp[lz] ) ;
131 break ;
132 }
133 default: {
134 cout << "Sym_tensor_trans::poisson() : problem with the grid!"
135 << endl ;
136 abort() ;
137 break ;
138 }
139 }
140 }
141 if (gri.get_type_r(nz-1) == UNSURR)
142 bornes[nz] = 1./(alp[nz-1] + bet[nz-1]) ;
143 else
144 bornes[nz] = alp[nz-1] + bet[nz-1] ;
145
146 const Mg3d& gr2 = *gri.get_non_axi() ;
147 Map_af mp2(gr2, bornes) ;
148 int np2 = ( np > 3 ? np : 4 ) ;
149
150 Sym_tensor sou_cart(mp2, CON, mp2.get_bvect_spher()) ;
151 for (int l=1; l<=3; l++)
152 for (int c=l; c<=3; c++) {
153 switch (this->operator()(l,c).get_etat() ) {
154 case ETATZERO: {
155 sou_cart.set(l,c).set_etat_zero() ;
156 break ;
157 }
158 case ETATUN: {
159 sou_cart.set(l,c).set_etat_one() ;
160 break ;
161 }
162 case ETATQCQ : {
163 sou_cart.set(l,c).allocate_all() ;
164 for (int lz=0; lz<nz; lz++)
165 for (int k=0; k<np2; k++)
166 for (int j=0; j<nt; j++)
167 for(int i=0; i<gr2.get_nr(lz); i++)
168 sou_cart.set(l,c).set_grid_point(lz, k, j, i)
169 = this->operator()(l,c).val_grid_point(lz, 0, j, i) ;
170 break ;
171 }
172 default: {
173 cout <<
174 "Sym_tensor_trans::poisson() : source in undefined state!"
175 << endl ;
176 abort() ;
177 break ;
178 }
179 }
180 sou_cart.set(l,c).set_dzpuis(this->operator()(l,c).get_dzpuis()) ;
181 }
182 sou_cart.std_spectral_base() ;
183 sou_cart.change_triad(mp2.get_bvect_cart()) ;
184 Sym_tensor res_cart(mp2, CON, mp2.get_bvect_cart()) ;
185 for (int i=1; i<=3; i++)
186 for(int j=i; j<=3; j++)
187 res_cart.set(i,j) = sou_cart(i,j).poisson() ;
188 res_cart.change_triad(mp2.get_bvect_spher()) ;
189 Scalar res_A(*mp) ; Scalar big_A = res_cart.compute_A() ;
190 Scalar res_B(*mp) ; Scalar big_B = res_cart.compute_tilde_B_tt() ;
191
192 switch (big_A.get_etat() ) {
193 case ETATZERO: {
194 res_A.set_etat_zero() ;
195 break ;
196 }
197 case ETATUN : {
198 res_A.set_etat_one() ;
199 break ;
200 }
201 case ETATQCQ : {
202 res_A.allocate_all() ;
203 for (int lz=0; lz<nz; lz++)
204 for (int k=0; k<np; k++)
205 for (int j=0; j<nt; j++)
206 for(int i=0; i<gri.get_nr(lz); i++)
207 res_A.set_grid_point(lz, k, j, i)
208 = big_A.val_grid_point(lz, k, j, i) ;
209 break ;
210 }
211 default: {
212 cout <<
213 "Sym_tensor_trans::poisson() : res_A in undefined state!"
214 << endl ;
215 abort() ;
216 break ;
217 }
218 }
219 res_A.set_spectral_base(big_A.get_spectral_base()) ;
220 int dzA = big_A.get_dzpuis() ;
221 res_A.set_dzpuis(dzA) ;
222
223 switch (big_B.get_etat() ) {
224 case ETATZERO: {
225 res_B.set_etat_zero() ;
226 break ;
227 }
228 case ETATUN : {
229 res_B.set_etat_one() ;
230 break ;
231 }
232 case ETATQCQ : {
233 res_B.allocate_all() ;
234 for (int lz=0; lz<nz; lz++)
235 for (int k=0; k<np; k++)
236 for (int j=0; j<nt; j++)
237 for(int i=0; i<gri.get_nr(lz); i++)
238 res_B.set_grid_point(lz, k, j, i)
239 = big_B.val_grid_point(lz, k, j, i) ;
240 break ;
241 }
242 default: {
243 cout <<
244 "Sym_tensor_trans::poisson() : res_B in undefined state!"
245 << endl ;
246 abort() ;
247 break ;
248 }
249 }
250 res_B.set_spectral_base(big_B.get_spectral_base()) ;
251 int dzB = big_B.get_dzpuis() ;
252 res_B.set_dzpuis(dzB) ;
253
254 resu.set_AtBtt_det_one(res_A, res_B, h_guess) ;
255
256 delete [] bornes ;
257 }
258 else {
259 assert (np >=4) ;
261 sou_cart.change_triad(mp->get_bvect_cart()) ;
262
263 Sym_tensor res_cart(*mp, CON, mp->get_bvect_cart()) ;
264 for (int i=1; i<=3; i++)
265 for(int j=i; j<=3; j++)
266 res_cart.set(i,j) = sou_cart(i,j).poisson() ;
267
268 res_cart.change_triad(*triad) ;
269
270 resu.set_AtBtt_det_one(res_cart.compute_A(), res_cart.compute_tilde_B_tt(), h_guess) ;
271
272 }
273#ifndef NDEBUG
274 Vector dive = resu.divergence(*met_div) ;
275 dive.dec_dzpuis(2) ;
276 maxabs(dive, "Sym_tensor_trans::poisson : divergence of the solution") ;
277#endif
278 return resu ;
279}
280
281
282}
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Affine radial mapping.
Definition map.h:2027
Multi-domain grid.
Definition grilles.h:273
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Transverse symmetric tensors of rank 2.
Definition sym_tensor.h:608
Sym_tensor_trans poisson(const Scalar *h_guess=0x0) const
Computes the solution of a tensorial transverse Poisson equation with *this as a source:
const Metric *const met_div
Metric with respect to which the divergence and the trace are defined.
Definition sym_tensor.h:614
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Tensor field of valence 1.
Definition vector.h:188
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
const Scalar & operator()(const Itbl &ind) const
Returns the value of a component (read-only version).
Definition tensor.C:798
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tensor.h:303
Tbl maxabs(const Tensor &aa, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maxima in each domain of the absolute values of the tensor components.
Lorene prototypes.
Definition app_hor.h:64