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