LORENE
vector.C
1 /*
2  * Methods of class Vector
3  *
4  * (see file vector.h for documentation)
5  *
6  */
7 
8 /*
9  * Copyright (c) 2003 Eric Gourgoulhon & 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 as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
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 
29 
30 
31 
32 /*
33  * $Id: vector.C,v 1.32 2024/01/26 17:44:25 g_servignat Exp $
34  * $Log: vector.C,v $
35  * Revision 1.32 2024/01/26 17:44:25 g_servignat
36  * Updated the Pseudopolytrope_1D class to be consistent with the paper (i.e. with a GPP in the middle)
37  *
38  * Revision 1.31 2016/12/05 16:18:18 j_novak
39  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
40  *
41  * Revision 1.30 2014/10/13 08:53:44 j_novak
42  * Lorene classes and functions now belong to the namespace Lorene.
43  *
44  * Revision 1.29 2014/10/06 15:13:20 j_novak
45  * Modified #include directives to use c++ syntax.
46  *
47  * Revision 1.28 2008/10/29 14:09:14 jl_cornou
48  * Spectral bases for pseudo vectors and curl added
49  *
50  * Revision 1.27 2008/08/27 08:52:23 jl_cornou
51  * Added fonctions for angular potential A
52  *
53  * Revision 1.26 2007/12/21 16:07:08 j_novak
54  * Methods to filter Tensor, Vector and Sym_tensor objects.
55  *
56  * Revision 1.25 2005/02/14 13:01:50 j_novak
57  * p_eta and p_mu are members of the class Vector. Most of associated functions
58  * have been moved from the class Vector_divfree to the class Vector.
59  *
60  * Revision 1.24 2005/01/25 15:37:35 j_novak
61  * Solved some dzpuis problem...
62  *
63  * Revision 1.23 2005/01/12 16:48:23 j_novak
64  * Better treatment of the case where all vector components are null in
65  * decompose_div .
66  *
67  * Revision 1.22 2004/10/12 09:58:25 j_novak
68  * Better memory management.
69  *
70  * Revision 1.21 2004/10/11 09:46:31 j_novak
71  * Speed improvements.
72  *
73  * Revision 1.20 2004/05/09 20:55:05 e_gourgoulhon
74  * Added method flux.
75  *
76  * Revision 1.19 2004/03/29 11:57:45 e_gourgoulhon
77  * Added methods ope_killing and ope_killing_conf.
78  *
79  * Revision 1.18 2004/02/26 22:48:50 e_gourgoulhon
80  * -- Method divergence: call to Tensor::divergence and cast of the
81  * result.
82  * -- Added method derive_lie.
83  *
84  * Revision 1.17 2004/02/24 09:46:20 j_novak
85  * Correction to cope with SGI compiler's warnings.
86  *
87  * Revision 1.16 2004/02/20 10:53:04 j_novak
88  * Added the matching of the potential adapted to the case where the
89  * vector is the source of a Poisson equation (dzpuis = 4).
90  *
91  * Revision 1.15 2004/01/30 10:30:17 j_novak
92  * Changed dzpuis handling in Vector::decompose_div (this may be temporary).
93  *
94  * Revision 1.14 2003/12/30 23:09:47 e_gourgoulhon
95  * Change in methods derive_cov() and divergence() to take into account
96  * the change of name: Metric::get_connect() --> Metric::connect().
97  *
98  * Revision 1.13 2003/12/19 15:18:16 j_novak
99  * Shadow variables hunt
100  *
101  * Revision 1.12 2003/10/29 11:04:34 e_gourgoulhon
102  * dec2_dpzuis() replaced by dec_dzpuis(2).
103  * inc2_dpzuis() replaced by inc_dzpuis(2).
104  *
105  * Revision 1.11 2003/10/22 14:24:19 j_novak
106  * *** empty log message ***
107  *
108  * Revision 1.9 2003/10/20 13:00:38 j_novak
109  * Memory error corrected
110  *
111  * Revision 1.8 2003/10/20 09:32:12 j_novak
112  * Members p_potential and p_div_free of the Helmholtz decomposition
113  * + the method decompose_div(Metric).
114  *
115  * Revision 1.7 2003/10/16 14:21:37 j_novak
116  * The calculation of the divergence of a Tensor is now possible.
117  *
118  * Revision 1.6 2003/10/13 13:52:40 j_novak
119  * Better managment of derived quantities.
120  *
121  * Revision 1.5 2003/10/06 13:58:48 j_novak
122  * The memory management has been improved.
123  * Implementation of the covariant derivative with respect to the exact Tensor
124  * type.
125  *
126  * Revision 1.4 2003/10/05 21:14:20 e_gourgoulhon
127  * Added method std_spectral_base().
128  *
129  * Revision 1.3 2003/10/03 14:10:32 e_gourgoulhon
130  * Added constructor from Tensor.
131  *
132  * Revision 1.2 2003/10/03 14:08:46 j_novak
133  * Removed old change_trid...
134  *
135  * Revision 1.1 2003/09/26 08:05:31 j_novak
136  * New class Vector.
137  *
138  *
139  * $Header: /cvsroot/Lorene/C++/Source/Tensor/vector.C,v 1.32 2024/01/26 17:44:25 g_servignat Exp $
140  *
141  */
142 
143 // Headers C
144 #include <cstdlib>
145 #include <cassert>
146 #include <cmath>
147 
148 // Headers Lorene
149 #include "metric.h"
150 #include "proto.h"
151 #include "matrice.h"
152 #include "nbr_spx.h"
153 
154 
155  //--------------//
156  // Constructors //
157  //--------------//
158 
159 // Standard constructor
160 // --------------------
161 namespace Lorene {
162 Vector::Vector(const Map& map, int tipe, const Base_vect& triad_i)
163  : Tensor(map, 1, tipe, triad_i) {
164 
165  set_der_0x0() ;
166 
167 }
168 
169 // Standard constructor with the triad passed as a pointer
170 // -------------------------------------------------------
171 Vector::Vector(const Map& map, int tipe, const Base_vect* triad_i)
172  : Tensor(map, 1, tipe, *triad_i) {
173 
174  set_der_0x0() ;
175 }
176 
177 // Copy constructor
178 // ----------------
179 Vector::Vector (const Vector& source) :
180  Tensor(source) {
181 
182  assert(valence == 1) ;
183  set_der_0x0() ;
184 
185 }
186 
187 
188 // Constructor from a {\tt Tensor}.
189 //--------------------------------
190 Vector::Vector(const Tensor& uu) : Tensor(uu) {
191 
192  assert(valence == 1) ;
193  set_der_0x0() ;
194 
195 }
196 
197 
198 // Constructor from a file
199 // -----------------------
200 Vector::Vector(const Map& mapping, const Base_vect& triad_i, FILE* fd) :
201  Tensor(mapping, triad_i, fd) {
202 
203  assert ( (valence == 1) && (n_comp == 3) ) ;
204  set_der_0x0() ;
205 
206 }
207 
208 
209  //--------------//
210  // Destructor //
211  //--------------//
212 
213 
215 
217 
218 }
219 
220 
221  //-------------------//
222  // Memory managment //
223  //-------------------//
224 
225 void Vector::del_deriv() const {
226 
227  for (int i=0; i<N_MET_MAX; i++)
228  del_derive_met(i) ;
229 
230  if (p_A != 0x0) delete p_A ;
231  if (p_eta != 0x0) delete p_eta ;
232  if (p_mu != 0x0) delete p_mu ;
233  set_der_0x0() ;
235 
236 }
237 
238 void Vector::set_der_0x0() const {
239 
240  for (int i=0; i<N_MET_MAX; i++)
241  set_der_met_0x0(i) ;
242  p_A = 0x0 ;
243  p_eta = 0x0 ;
244  p_mu = 0x0 ;
245 
246 }
247 
248 void Vector::del_derive_met(int j) const {
249 
250  assert( (j>=0) && (j<N_MET_MAX) ) ;
251 
252  if (met_depend[j] != 0x0) {
253  if (p_potential[j] != 0x0)
254  delete p_potential[j] ;
255  if (p_div_free[j] != 0x0)
256  delete p_div_free[j] ;
257 
258  set_der_met_0x0(j) ;
259 
261  }
262 }
263 
264 void Vector::set_der_met_0x0(int i) const {
265 
266  assert( (i>=0) && (i<N_MET_MAX) ) ;
267 
268  p_potential[i] = 0x0 ;
269  p_div_free[i] = 0x0 ;
270 
271 }
272 
273 void Vector::operator=(const Vector& t) {
274 
275  triad = t.triad ;
276 
277  assert(t.type_indice(0) == type_indice(0)) ;
278 
279  for (int i=0 ; i<3 ; i++) {
280  *cmp[i] = *t.cmp[i] ;
281  }
282  del_deriv() ;
283 }
284 
285 void Vector::operator=(const Tensor& t) {
286 
287  assert (t.valence == 1) ;
288 
289  triad = t.triad ;
290 
291  assert(t.type_indice(0) == type_indice(0)) ;
292 
293  for (int i=0 ; i<3 ; i++) {
294  *cmp[i] = *t.cmp[i] ;
295  }
296  del_deriv() ;
297 }
298 
299 
300 
301 // Affectation d'une composante :
302 Scalar& Vector::set(int index) {
303 
304  assert ( (index>=1) && (index<=3) ) ;
305 
306  del_deriv() ;
307 
308  return *cmp[index - 1] ;
309 }
310 
311 const Scalar& Vector::operator()(int index) const {
312 
313  assert ( (index>=1) && (index<=3) ) ;
314 
315  return *cmp[index - 1] ;
316 
317 }
318 
319 
320 // Sets the standard spectal bases of decomposition for each component
321 
323 
324  Base_val** bases = 0x0 ;
325 
326  if ( triad->identify() == (mp->get_bvect_cart()).identify() ) {
327 
328  // Cartesian case
329  bases = mp->get_mg()->std_base_vect_cart() ;
330 
331  }
332  else {
333  // Spherical case
334  assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
335  bases = mp->get_mg()->std_base_vect_spher() ;
336  }
337 
338  for (int i=0 ; i<3 ; i++) {
339  cmp[i]->set_spectral_base( *bases[i] ) ;
340  }
341 
342  for (int i=0 ; i<3 ; i++) {
343  delete bases[i] ;
344  }
345  delete [] bases ;
346 
347 
348 }
349 
350 // Sets the standard spectral bases of decomposition for each component for a pseudo vector
351 
353 
354  Base_val** bases = 0x0 ;
355 
356  if ( triad->identify() == (mp->get_bvect_cart()).identify() ) {
357 
358  // Cartesian case
359  bases = mp->get_mg()->pseudo_base_vect_cart() ;
360 
361  }
362  else {
363  // Spherical case
364  assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
365  bases = mp->get_mg()->pseudo_base_vect_spher() ;
366  }
367 
368  for (int i=0 ; i<3 ; i++) {
369  cmp[i]->set_spectral_base( *bases[i] ) ;
370  }
371 
372  for (int i=0 ; i<3 ; i++) {
373  delete bases[i] ;
374  }
375  delete [] bases ;
376 
377 
378 }
379 
380 
381 
382  //-------------------------------//
383  // Computational methods //
384  //-------------------------------//
385 
386 
387 const Scalar& Vector::divergence(const Metric& metre) const {
388 
389  const Scalar* pscal =
390  dynamic_cast<const Scalar*>( &(Tensor::divergence(metre)) ) ;
391 
392  assert(pscal != 0x0) ;
393 
394  return *pscal ;
395 }
396 
397 
398 Vector Vector::derive_lie(const Vector& vv) const {
399 
400  Vector resu(*mp, type_indice(0), triad) ;
401 
402  compute_derive_lie(vv, resu) ;
403 
404  return resu ;
405 
406 }
407 
409 
410  if ( triad->identify() == (mp->get_bvect_cart()).identify() ) {
411  const Metric_flat& metc = mp->flat_met_cart() ;
412  Vector_divfree resu(*mp, mp->get_bvect_cart(), metc) ;
413  resu.set(1)= cmp[2]->dsdy() - cmp[1]->dsdz();
414  resu.set(2)= cmp[0]->dsdz() - cmp[2]->dsdx();
415  resu.set(3)= cmp[1]->dsdx() - cmp[0]->dsdy();
416  resu.pseudo_spectral_base();
417  return resu ;
418  }
419  else {
420  assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
421  const Metric_flat& mets = mp->flat_met_spher() ;
422  Vector_divfree resu(*mp, mp->get_bvect_spher(), mets);
423  Scalar tmp = *cmp[2] ;
424  tmp.div_tant();
425  tmp += cmp[2]->dsdt();
426  tmp.div_r();
427  resu.set(1) = tmp - cmp[1]->srstdsdp() ;
428  tmp = *cmp[2] ;
429  tmp.mult_r();
430  tmp = tmp.dsdr();
431  tmp.div_r();
432  resu.set(2) = cmp[0]->srstdsdp() - tmp ;
433  tmp = *cmp[1];
434  tmp.mult_r();
435  resu.set(3) = tmp.dsdr() - cmp[0]->dsdt() ;
436  resu.set(3).div_r();
437  resu.pseudo_spectral_base();
438  return resu ;
439  }
440 
441 }
442 
443 
445 
446  Sym_tensor resu(*mp, type_indice(0), *triad) ;
447 
448  if (type_indice(0) == CON ) {
449  for (int i=1; i<=3; i++) {
450  for (int j=i; j<=3; j++) {
451  resu.set(i,j) = derive_con(gam)(i,j) + derive_con(gam)(j,i) ;
452  }
453  }
454  }
455  else {
456  for (int i=1; i<=3; i++) {
457  for (int j=i; j<=3; j++) {
458  resu.set(i,j) = derive_cov(gam)(i,j) + derive_cov(gam)(j,i) ;
459  }
460  }
461  }
462 
463  return resu ;
464 
465 }
466 
467 
469 
470  Sym_tensor resu(*mp, type_indice(0), *triad) ;
471 
472  if (type_indice(0) == CON ) {
473  for (int i=1; i<=3; i++) {
474  for (int j=i; j<=3; j++) {
475  resu.set(i,j) = derive_con(gam)(i,j) + derive_con(gam)(j,i)
476  - 0.6666666666666666* divergence(gam)
477  * gam.con()(i,j) ;
478  }
479  }
480  }
481  else {
482  for (int i=1; i<=3; i++) {
483  for (int j=i; j<=3; j++) {
484  resu.set(i,j) = derive_cov(gam)(i,j) + derive_cov(gam)(j,i)
485  - 0.6666666666666666* derive_con(gam).trace()
486  * gam.cov()(i,j) ;
487  }
488  }
489  }
490 
491  return resu ;
492 
493 }
494 
495 
496 
497 
498 const Scalar& Vector::potential(const Metric& metre) const {
499 
500  set_dependance(metre) ;
501  int j = get_place_met(metre) ;
502  assert ((j>=0) && (j<N_MET_MAX)) ;
503  if (p_potential[j] == 0x0) {
504  decompose_div(metre) ;
505  }
506 
507  return *p_potential[j] ;
508 }
509 
510 const Vector_divfree& Vector::div_free(const Metric& metre) const {
511 
512  set_dependance(metre) ;
513  int j = get_place_met(metre) ;
514  assert ((j>=0) && (j<N_MET_MAX)) ;
515  if (p_div_free[j] == 0x0) {
516  decompose_div(metre) ;
517  }
518 
519  return *p_div_free[j] ;
520 }
521 
522 void Vector::decompose_div(const Metric& metre) const {
523 
524  assert( type_indice(0) == CON ) ; //Only for contravariant vectors...
525 
526  set_dependance(metre) ;
527  int ind = get_place_met(metre) ;
528  assert ((ind>=0) && (ind<N_MET_MAX)) ;
529 
530  if ( (p_potential[ind] != 0x0) && (p_div_free[ind] != 0x0) )
531  return ; // Nothing to do ...
532 
533  else {
534  if (p_div_free[ind] != 0x0)
535  delete p_div_free[ind] ;
536  if (p_potential[ind] != 0x0)
537  delete p_potential[ind] ;
538 
539  const Mg3d* mmg = mp->get_mg() ;
540  int nz = mmg->get_nzone() ;
541 
542  int dzp = cmp[0]->get_dzpuis() ;
543 #ifndef NDEBUG
544  bool dz_zero = cmp[0]->check_dzpuis(0) ;
545  bool dz_four = cmp[0]->check_dzpuis(4) ;
546 #endif
547  assert( dz_zero || dz_four) ;
548  assert (cmp[1]->check_dzpuis(dzp)) ;
549  assert (cmp[2]->check_dzpuis(dzp)) ;
550 
551  Scalar dive = divergence(metre) ;
552 
553  if (dive.get_etat() == ETATZERO) {
554  p_potential[ind] = new Scalar(*mp) ;
555  p_potential[ind]->set_etat_zero() ;
556  p_potential[ind]->set_dzpuis(dzp) ;
557  }
558  else {
559  //No matching of the solution at this point
560  p_potential[ind] = new Scalar(dive.poisson()) ;
561 
562  if (dzp == 4) {
563  assert (mmg->get_type_r(nz-1) == UNSURR) ;
564  // Let's now do the matching ... in the case dzpuis = 4
565  const Map_af* mapping = dynamic_cast<const Map_af*>(mp) ;
566  assert (mapping != 0x0) ;
567  Valeur& val = p_potential[ind]->set_spectral_va() ;
568  val.ylm() ;
569  Mtbl_cf& mtc = *val.c_cf ;
570  if (nz > 1) {
571  int np = mmg->get_np(0) ;
572  int nt = mmg->get_nt(0) ;
573 #ifndef NDEBUG
574  for (int lz=0; lz<nz; lz++) {
575  assert (mmg->get_np(lz) == np) ;
576  assert (mmg->get_nt(lz) == nt) ;
577  }
578 #endif
579  int lmax = 0 ;
580  for (int k=0; k<np+1; k++)
581  for (int j=0; j<nt; j++)
582  if ( nullite_plm(j, nt, k, np, val.base)) {
583  int m_quant, l_quant, base_r ;
584  donne_lm (nz, 0, j, k, val.base, m_quant,
585  l_quant, base_r) ;
586  lmax = (l_quant > lmax ? l_quant : lmax) ;
587  }
588  Scalar** ri = new Scalar*[lmax+1] ;
589  Scalar** rmi = new Scalar*[lmax+1] ;
590  Scalar erre(*mp) ;
591  erre = mp->r ;
592  for (int l=0; l<=lmax; l++) {
593  ri[l] = new Scalar(*mp) ;
594  rmi[l] = new Scalar(*mp) ;
595  if (l == 0) *(ri[l]) = 1. ;
596  else *(ri[l]) = pow(erre, l) ;
597  ri[l]->annule_domain(nz - 1) ;
598  ri[l]->std_spectral_base() ; //exact base in r will be set later
599  if (l==2) *(rmi[l]) = 1 ;
600  else *(rmi[l]) = pow(erre, 2-l) ;
601  rmi[l]->annule(0,nz-2) ;
602  Scalar tmp = pow(erre, -l-1) ;
603  tmp.annule_domain(nz-1) ;
604  tmp.annule_domain(0) ;
605  *(rmi[l]) += tmp ;
606  rmi[l]->set_dzpuis(3) ;
607  rmi[l]->std_spectral_base() ;//exact base in r will be set later
608  }
609 
610  for (int k=0; k<np+1; k++) {
611  for (int j=0; j<nt; j++) {
612  if ( nullite_plm(j, nt, k, np, val.base)) {
613  int m_quant, l_quant, base_r, l, c ;
614  donne_lm (nz, 0, j, k, val.base, m_quant, l_quant, base_r) ;
615  bool lzero = (l_quant == 0) || (l_quant == 1) ;
616  int nb_hom_ced = (l_quant < 2 ? 0 : 1) ;
617 
618  int taille = 2*(nz-1) - 1 + nb_hom_ced ;
619  Tbl deuz(taille) ;
620  deuz.set_etat_qcq() ;
621  Matrice systeme(taille,taille) ;
622  systeme.set_etat_qcq() ;
623  for (l=0; l<taille; l++)
624  for (c=0; c<taille; c++) systeme.set(l,c) = 0 ;
625  for (l=0; l<taille; l++) deuz.set(l) = 0 ;
626 
627  //---------
628  // Nucleus
629  //---------
630  assert(mmg->get_type_r(0) == RARE) ;
631  assert ((base_r == R_CHEBP)||(base_r == R_CHEBI)) ;
632  ri[l_quant]->set_spectral_va().set_base_r(0, base_r) ;
633 
634  double alpha = mapping->get_alpha()[0] ;
635  int nr = mmg->get_nr(0) ;
636  Tbl partn(nr) ;
637  partn.set_etat_qcq() ;
638  for (int i=0; i<nr; i++)
639  partn.set(i) = mtc(0, k, j, i) ;
640  l=0 ; c=0 ;
641 
642  systeme.set(l,c) += pow(alpha, double(l_quant)) ;
643 
644  deuz.set(l) -= val1_dern_1d(0, partn, base_r) ;
645  l++ ;
646 
647  if (!lzero) {
648  systeme.set(l,c) += l_quant * pow(alpha, double(l_quant - 1)) ;
649 
650  deuz.set(l) -= val1_dern_1d(1, partn, base_r) / alpha ;
651  }
652 
653  //----------
654  // Shells
655  //----------
656 
657  for (int lz=1; lz<nz-1; lz++) {
658  alpha = mapping->get_alpha()[lz] ;
659  double beta = mapping->get_beta()[lz] ;
660  double xm1 = beta - alpha ;
661  double xp1 = alpha + beta ;
662  nr = mmg->get_nr(lz) ;
663  Tbl parts(nr) ;
664  parts.set_etat_qcq() ;
665  for (int i=0; i<nr; i++)
666  parts.set(i) = mtc(lz, k, j, i) ;
667 
668  //Function at x = -1
669  l-- ;
670  c = 2*lz - 1 ;
671  systeme.set(l,c) -= pow(xm1, double(l_quant)) ;
672  c++;
673  systeme.set(l,c) -= pow(xm1, double(-l_quant-1)) ;
674 
675  deuz.set(l) += valm1_dern_1d(0, parts, R_CHEB) ;
676 
677  if ((lz>1) || (!lzero)) {
678  //First derivative at x=-1
679  l++ ;
680  c-- ;
681  systeme.set(l,c) -= l_quant*pow(xm1, double(l_quant - 1)) ;
682  c++;
683  systeme.set(l,c) -= (-l_quant - 1)*
684  pow(xm1, double(-l_quant - 2)) ;
685 
686  deuz.set(l) += valm1_dern_1d(1, parts, R_CHEB) / alpha ;
687  }
688 
689  //Function at x = 1
690  l++ ;
691  c-- ;
692  systeme.set(l,c) += pow(xp1, double(l_quant)) ;
693  c++;
694  systeme.set(l,c) += pow(xp1, double(-l_quant-1)) ;
695 
696  deuz.set(l) -= val1_dern_1d(0, parts, R_CHEB) ;
697 
698  //First derivative at x = 1
699  l++ ;
700  c-- ;
701  systeme.set(l,c) += l_quant*pow(xp1, double(l_quant - 1)) ;
702  c++;
703  systeme.set(l,c) += (-l_quant - 1)*
704  pow(xp1, double(-l_quant - 2)) ;
705 
706  deuz.set(l) -= val1_dern_1d(1, parts, R_CHEB) / alpha ;
707 
708  } //End of the loop on different shells
709 
710  //-------------------------------
711  // Compactified external domain
712  //-------------------------------
713  assert(mmg->get_type_r(nz-1) == UNSURR) ;
714  nr = mmg->get_nr(nz-1) ;
715  Tbl partc(nr) ;
716  partc.set_etat_qcq() ;
717  for (int i=0; i<nr; i++)
718  partc.set(i) = mtc(nz-1, k, j, i) ;
719 
720  alpha = mapping->get_alpha()[nz-1] ;
721  double beta = mapping->get_beta()[nz-1] ;
722  double xm1 = beta - alpha ; // 1 / r_left
723  double part0, part1 ;
724  part0 = valm1_dern_1d(0, partc, R_CHEB) ;
725  part1 = xm1*valm1_dern_1d(1, partc, R_CHEB) / alpha ;
726  assert (p_potential[ind]->get_dzpuis() == 3) ;
727 
728  //Function at x = -1
729  l--;
730  if (!lzero) {
731  c++;
732  systeme.set(l,c) -= pow(xm1, double(l_quant+1)) ;
733  }
734  deuz.set(l) += part0*xm1*xm1*xm1 ;
735 
736  // First derivative at x = -1
737  l++ ;
738  if (!lzero) {
739  systeme.set(l,c) -= (-l_quant - 1)*
740  pow(xm1, double(l_quant + 2)) ;
741  }
742  if ( (nz > 2) || (!lzero))
743  deuz.set(l) += -xm1*xm1*xm1*xm1*(3*part0 + part1) ;
744 
745  //--------------------------------------
746  // Solution of the linear system
747  //--------------------------------------
748 
749  int inf = 1 + nb_hom_ced;
750  int sup = 3 - nb_hom_ced;
751  systeme.set_band(sup, inf) ;
752  systeme.set_lu() ;
753  Tbl facteur(systeme.inverse(deuz)) ;
754  ri[l_quant]->set_spectral_va().coef() ;
755  rmi[l_quant]->set_spectral_va().coef() ;
756 
757  //Linear combination in the nucleus
758  nr = mmg->get_nr(0) ;
759  for (int i=0; i<nr; i++)
760  mtc.set(0, k, j, i) +=
761  facteur(0)*(*(ri[l_quant]->get_spectral_va().c_cf))(0, 0, 0, i) ;
762 
763  //Linear combination in the shells
764  for (int lz=1; lz<nz-1; lz++) {
765  nr = mmg->get_nr(lz) ;
766  for (int i=0; i<nr; i++)
767  mtc.set(lz, k, j, i) += facteur(2*lz - 1)*
768  (*(ri[l_quant]->get_spectral_va().c_cf))(lz, 0, 0, i) ;
769  for (int i=0; i<nr; i++)
770  mtc.set(lz, k, j, i) += facteur(2*lz)*
771  (*(rmi[l_quant]->get_spectral_va().c_cf))(lz, 0, 0, i) ;
772  }
773 
774  //Linear combination in the CED
775  nr = mmg->get_nr(nz-1) ;
776  if (!lzero) {
777  for (int i=0; i<nr; i++)
778  mtc.set(nz-1, k, j, i) +=
779  facteur(taille - 1)*
780  (*(rmi[l_quant]->get_spectral_va().c_cf))(nz-1, 0, 0, i) ;
781  }
782  } //End of nullite_plm ...
783 
784  } //End of j/theta loop
785  } //End of k/phi loop
786 
787  for (int l=0; l<=lmax; l++) {
788  delete ri[l] ;
789  delete rmi[l] ;
790  }
791  delete [] ri ;
792  delete [] rmi ;
793 
794  } //End of the case of more than one domain
795 
796  val.ylm_i() ;
797 
798  } //End of the case dzp = 4
799  }
800  p_div_free[ind] = new Vector_divfree(*mp, *triad, metre) ;
801 
802  Vector gradient = p_potential[ind]->derive_con(metre) ;
803  if (dzp != 4) gradient.dec_dzpuis(2) ;
804 
805  *p_div_free[ind] = ( *this - gradient ) ;
806 
807  }
808 
809 }
810 
811 
812 
813 double Vector::flux(double radius, const Metric& met) const {
814 
815  assert(type_indice(0) == CON) ;
816 
817  const Map_af* mp_af = dynamic_cast<const Map_af*>(mp) ;
818  if (mp_af == 0x0) {
819  cerr <<
820  "Vector::flux : the case of a mapping outside the class Map_af\n"
821  << " is not implemented yet !" << endl ;
822  abort() ;
823  }
824 
825  const Metric_flat* ff = dynamic_cast<const Metric_flat*>(&met) ;
826  if (ff == 0x0) {
827  cerr <<
828  "Vector::flux : the case of a non flat metric is not implemented yet !"
829  << endl ;
830  abort() ;
831  }
832 
833  const Base_vect_cart* bcart = dynamic_cast<const Base_vect_cart*>(triad) ;
834  Vector* vspher = 0x0 ;
835  if (bcart != 0x0) { // switch to spherical components:
836  vspher = new Vector(*this) ;
837  vspher->change_triad(mp->get_bvect_spher()) ;
838  }
839  else assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ;
840 
841  const Vector* vv = (bcart != 0x0) ? vspher : this ;
842 
843  const Scalar& vr = vv->operator()(1) ;
844 
845  double resu ;
846  if (radius == __infinity) {
847  resu = mp_af->integrale_surface_infini(vr) ;
848  }
849  else {
850  resu = mp_af->integrale_surface(vr, radius) ;
851  }
852 
853  return resu ;
854 }
855 
856 void Vector::exponential_filter_r(int lzmin, int lzmax, int p,
857  double alpha) {
858  if( triad->identify() == (mp->get_bvect_cart()).identify() )
859  for (int i=0; i<n_comp; i++)
860  cmp[i]->exponential_filter_r(lzmin, lzmax, p, alpha) ;
861  else {
862  assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
863  Scalar vr_tmp = operator()(1) ;
864  vr_tmp.exponential_filter_r(lzmin, lzmax, p, alpha) ;
865  Scalar eta_tmp = eta() ;
866  eta_tmp.exponential_filter_r(lzmin, lzmax, p, alpha) ;
867  Scalar mu_tmp = mu() ;
868  mu_tmp.exponential_filter_r(lzmin, lzmax, p, alpha) ;
869  set_vr_eta_mu(vr_tmp, eta_tmp, mu_tmp) ;
870  }
871 }
872 
873 void Vector::exponential_filter_ylm(int lzmin, int lzmax, int p,
874  double alpha) {
875  if( triad->identify() == (mp->get_bvect_cart()).identify() )
876  for (int i=0; i<n_comp; i++)
877  cmp[i]->exponential_filter_ylm(lzmin, lzmax, p, alpha) ;
878  else {
879  assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
880  Scalar vr_tmp = operator()(1) ;
881  vr_tmp.exponential_filter_ylm(lzmin, lzmax, p, alpha) ;
882  Scalar eta_tmp = eta() ;
883  eta_tmp.exponential_filter_ylm(lzmin, lzmax, p, alpha) ;
884  Scalar mu_tmp = mu() ;
885  mu_tmp.exponential_filter_ylm(lzmin, lzmax, p, alpha) ;
886  set_vr_eta_mu(vr_tmp, eta_tmp, mu_tmp) ;
887  }
888 }
889 }
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
Metric for tensor calculation.
Definition: metric.h:90
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:293
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
double flux(double radius, const Metric &met) const
Computes the flux of the vector accross a sphere r = const.
Definition: vector.C:813
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
Scalar * p_potential[N_MET_MAX]
The potential giving the gradient part in the Helmholtz decomposition of any 3D vector ...
Definition: vector.h:198
Sym_tensor ope_killing(const Metric &gam) const
Computes the Killing operator associated with a given metric.
Definition: vector.C:444
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:134
int n_comp
Number of stored components, depending on the symmetry.
Definition: tensor.h:318
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:151
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
const Vector_divfree & div_free(const Metric &) const
Returns the div-free vector in the Helmholtz decomposition.
Definition: vector.C:510
double integrale_surface_infini(const Cmp &ci) const
Performs the surface integration of ci at infinity.
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:397
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition: matrice.C:395
virtual void del_derive_met(int) const
Logical destructor of the derivatives depending on the i-th element of met_depend ...
Definition: tensor.C:423
Lorene prototypes.
Definition: app_hor.h:67
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:304
void set_dependance(const Metric &) const
To be used to describe the fact that the derivatives members have been calculated with met ...
Definition: tensor.C:462
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition: matrice.C:427
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:141
const Scalar & dsdt() const
Returns of *this .
Definition: scalar_deriv.C:208
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Flat metric for tensor calculation.
Definition: metric.h:261
const Scalar & dsdz() const
Returns of *this , where .
Definition: scalar_deriv.C:328
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: scalar_pde.C:139
Vector_divfree * p_div_free[N_MET_MAX]
The divergence-free vector of the Helmholtz decomposition of any 3D vector .
Definition: vector.h:205
Scalar * p_A
Field defined by Insensitive to the longitudinal part of the vector, related to the curl...
Definition: vector.h:241
Base class for coordinate mappings.
Definition: map.h:688
const Scalar & dsdx() const
Returns of *this , where .
Definition: scalar_deriv.C:266
virtual void exponential_filter_r(int lzmin, int lzmax, int p, double alpha=-16.)
Applies an exponential filter to the spectral coefficients in the radial direction.
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tensor.h:309
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
void set_der_0x0() const
Sets the pointers on derived quantities to 0x0.
Definition: vector.C:238
int get_place_met(const Metric &) const
Returns the position of the pointer on metre in the array met_depend .
Definition: tensor.C:452
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition: scalar_deriv.C:402
Tensor field of valence 1.
Definition: vector.h:188
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition: base_vect.h:105
Vector(const Map &map, int tipe, const Base_vect &triad_i)
Standard constructor.
Definition: vector.C:162
const Metric_flat & flat_met_cart() const
Returns the flat metric associated with the Cartesian coordinates and with components expressed in th...
Definition: map.C:334
Base_val ** pseudo_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a pseudo-vector.
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:814
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:322
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition: tensor.C:817
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
Sym_tensor ope_killing_conf(const Metric &gam) const
Computes the conformal Killing operator associated with a given metric.
Definition: vector.C:468
virtual void exponential_filter_r(int lzmin, int lzmax, int p, double alpha=-16.)
Applies exponential filters to all components (see Scalar::exponential_filter_r ).
Definition: vector.C:856
Base_val ** pseudo_base_vect_spher() const
Returns the standard spectral bases for the spherical components of a pseudo-vector.
void div_r()
Division by r everywhere; dzpuis is not changed.
Itbl type_indice
1D array of integers (class Itbl ) of size valence containing the type of each index: COV for a cova...
Definition: tensor.h:316
virtual void exponential_filter_ylm(int lzmin, int lzmax, int p, double alpha=-16.)
Applies exponential filters to all components (see Scalar::exponential_filter_ylm )...
Definition: vector.C:873
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
virtual void del_derive_met(int) const
Logical destructor of the derivatives depending on the i-th element of met_depend in the class Vector...
Definition: vector.C:248
void compute_derive_lie(const Vector &v, Tensor &resu) const
Computes the Lie derivative of this with respect to some vector field v (protected method; the public...
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written: .
Definition: vector_etamu.C:72
const Scalar & dsdy() const
Returns of *this , where .
Definition: scalar_deriv.C:297
Matrix handling.
Definition: matrice.h:152
virtual void del_deriv() const
Deletes the derived quantities.
Definition: vector.C:225
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:315
void set_vr_eta_mu(const Scalar &vr_i, const Scalar &eta_i, const Scalar &mu_i)
Defines the components through potentials and (see members p_eta and p_mu ), as well as the compon...
Definition: vector_etamu.C:198
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:563
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition: tensor.h:321
Base_val ** std_base_vect_spher() const
Returns the standard spectral bases for the spherical components of a vector.
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:387
virtual int identify() const =0
Returns a number to identify the sub-classe of Base_vect the object belongs to.
Scalar * p_mu
Field such that the angular components of the vector are written: .
Definition: vector.h:233
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Tensor handling.
Definition: tensor.h:294
void decompose_div(const Metric &) const
Makes the Helmholtz decomposition (see documentation of p_potential ) of this with respect to a given...
Definition: vector.C:522
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor.C:1011
double & set(int j, int i)
Read/write of a particuliar element.
Definition: matrice.h:277
const Tensor & divergence(const Metric &gam) const
Computes the divergence of this with respect to some metric .
Definition: tensor.C:1064
virtual void del_deriv() const
Deletes the derived quantities.
Definition: tensor.C:407
Cartesian vectorial bases (triads).
Definition: base_vect.h:201
void set_band(int up, int low) const
Calculate the band storage of *std.
Definition: matrice.C:367
void set_der_met_0x0(int) const
Sets all the i-th components of met_depend in the class Vector (p_potential , etc...) to 0x0.
Definition: vector.C:264
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
Definition: scalar.C:803
int valence
Valence of the tensor (0 = scalar, 1 = vector, etc...)
Definition: tensor.h:304
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:283
virtual void operator=(const Vector &a)
Assignment from a Vector.
Definition: vector.C:273
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition: tensor.C:1023
Multi-domain grid.
Definition: grilles.h:279
Vector derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition: vector.C:398
Bases of the spectral expansions.
Definition: base_val.h:325
void set_base_r(int l, int base_r)
Sets the expansion basis for r ( ) functions in a given domain.
Definition: valeur.C:839
Affine radial mapping.
Definition: map.h:2048
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
const Metric * met_depend[N_MET_MAX]
Array on the Metric &#39;s which were used to compute derived quantities, like p_derive_cov ...
Definition: tensor.h:333
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
const Vector_divfree curl() const
The curl of this with respect to a (flat) Metric .
Definition: vector.C:408
const Scalar & potential(const Metric &) const
Returns the potential in the Helmholtz decomposition.
Definition: vector.C:498
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:196
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written: .
Definition: vector_etamu.C:107
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:178
Basic array class.
Definition: tbl.h:164
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
Scalar * p_eta
Field such that the angular components of the vector are written: .
Definition: vector.h:219
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
Divergence-free vectors.
Definition: vector.h:724
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:491
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition: scalar.C:879
virtual void exponential_filter_ylm(int lzmin, int lzmax, int p, double alpha=-16.)
Applies an exponential filter to the spectral coefficients in the angular directions.
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:301
const Scalar & operator()(int) const
Readonly access to a component.
Definition: vector.C:311
void div_tant()
Division by .
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
Definition: map.C:324
virtual ~Vector()
Destructor.
Definition: vector.C:214
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607
const Scalar & srstdsdp() const
Returns of *this .
Definition: scalar_deriv.C:177
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
Coord r
r coordinate centered on the grid
Definition: map.h:736
virtual void pseudo_spectral_base()
Sets the standard spectal bases of decomposition for each component for a pseudo_vector.
Definition: vector.C:352
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166