LORENE
tensor.C
1/*
2 * Methods of class Tensor
3 *
4 * (see file tensor.h for documentation)
5 *
6 */
7
8/*
9 * Copyright (c) 2003-2004 Eric Gourgoulhon & Jerome Novak
10 *
11 * Copyright (c) 1999-2001 Philippe Grandclement (for preceding class Tenseur)
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 tensor_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/tensor.C,v 1.43 2014/10/13 08:53:44 j_novak Exp $" ;
33
34/*
35 * $Id: tensor.C,v 1.43 2014/10/13 08:53:44 j_novak Exp $
36 * $Log: tensor.C,v $
37 * Revision 1.43 2014/10/13 08:53:44 j_novak
38 * Lorene classes and functions now belong to the namespace Lorene.
39 *
40 * Revision 1.42 2014/10/06 15:13:20 j_novak
41 * Modified #include directives to use c++ syntax.
42 *
43 * Revision 1.41 2013/06/05 15:43:49 j_novak
44 * Suppression of leg_spectral_base()
45 *
46 * Revision 1.40 2013/01/11 15:44:54 j_novak
47 * Addition of Legendre bases (part 2).
48 *
49 * Revision 1.39 2007/12/21 16:07:08 j_novak
50 * Methods to filter Tensor, Vector and Sym_tensor objects.
51 *
52 * Revision 1.38 2005/10/25 08:56:38 p_grandclement
53 * addition of std_spectral_base in the case of odd functions near the origin
54 *
55 * Revision 1.37 2005/05/18 11:45:44 j_novak
56 * Added del_deriv() calls at the end of inc/dec_dzpuis.
57 *
58 * Revision 1.36 2004/07/08 12:21:53 j_novak
59 * Replaced tensor::annule_extern_c2 with tensor::annule_extern_cn for a
60 * more general transition.
61 *
62 * Revision 1.35 2004/06/17 06:55:07 e_gourgoulhon
63 * Added method annule_extern_c2.
64 *
65 * Revision 1.34 2004/03/04 09:47:51 e_gourgoulhon
66 * Method annule_domain(int, int): added call to virtual function
67 * del_deriv() at the end !
68 *
69 * Revision 1.33 2004/02/27 21:14:27 e_gourgoulhon
70 * Modif of method derive_con to create proper type of the result.
71 *
72 * Revision 1.32 2004/02/19 22:10:14 e_gourgoulhon
73 * Added argument "comment" in method spectral_display.
74 *
75 * Revision 1.31 2004/02/05 15:03:47 e_gourgoulhon
76 * Corrected bug in method derive_con().
77 *
78 * Revision 1.30 2004/01/27 13:05:11 j_novak
79 * Removed the method Tensor::mult_r_ced()
80 *
81 * Revision 1.29 2004/01/19 16:32:13 e_gourgoulhon
82 * Added operator()(int, int, int, int) and set(int, int, int, int)
83 * for direct access to components of valence 4 tensors.
84 *
85 * Revision 1.28 2004/01/04 20:55:23 e_gourgoulhon
86 * Method spectral_display(): added printing of type of class (through typeid).
87 *
88 * Revision 1.27 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.26 2003/12/27 15:01:50 e_gourgoulhon
93 * Method derive_con(): the position of the derivation index has
94 * been changed from the first one to the last one.
95 *
96 * Revision 1.25 2003/11/03 22:34:41 e_gourgoulhon
97 * Method dec_dzpuis: changed the name of argument dec --> decrem
98 * (in order not to shadow some globally defined dec).
99 *
100 * Revision 1.24 2003/10/29 11:02:54 e_gourgoulhon
101 * Functions dec_dzpuis and inc_dzpuis have now an integer argument to
102 * specify by which amount dzpuis is to be increased.
103 * Accordingly methods dec2_dzpuis and inc2_dzpuis have been suppressed
104 *
105 * Revision 1.23 2003/10/27 10:49:48 e_gourgoulhon
106 * Slightly modified operator<<.
107 *
108 * Revision 1.22 2003/10/19 19:57:00 e_gourgoulhon
109 * -- Added new method spectral_display
110 * -- slightly rearranged the operator<<
111 *
112 * Revision 1.21 2003/10/16 15:27:46 e_gourgoulhon
113 * Name of method annule(int ) changed to annule_domain(int ).
114 *
115 * Revision 1.20 2003/10/16 14:21:36 j_novak
116 * The calculation of the divergence of a Tensor is now possible.
117 *
118 * Revision 1.19 2003/10/13 13:52:40 j_novak
119 * Better managment of derived quantities.
120 *
121 * Revision 1.18 2003/10/11 16:47:10 e_gourgoulhon
122 * Suppressed the call to Ibtl::set_etat_qcq() after the construction
123 * of the Itbl's, thanks to the new property of the Itbl class.
124 *
125 * Revision 1.17 2003/10/11 14:48:40 e_gourgoulhon
126 * Line 344: suppressed assert(resu == -1)
127 * and added a break command to quit the loop.
128 *
129 * Revision 1.16 2003/10/08 14:24:09 j_novak
130 * replaced mult_r_zec with mult_r_ced
131 *
132 * Revision 1.15 2003/10/07 15:25:38 e_gourgoulhon
133 * Added call to del_derive_met in del_deriv().
134 *
135 * Revision 1.14 2003/10/07 09:10:00 j_novak
136 * Use of ::contract instead of up()
137 *
138 * Revision 1.13 2003/10/06 20:51:43 e_gourgoulhon
139 * In methods set: changed name "indices" to "idx" to avoid shadowing
140 * of class member.
141 *
142 * Revision 1.12 2003/10/06 16:17:31 j_novak
143 * Calculation of contravariant derivative and Ricci scalar.
144 *
145 * Revision 1.11 2003/10/06 13:58:48 j_novak
146 * The memory management has been improved.
147 * Implementation of the covariant derivative with respect to the exact Tensor
148 * type.
149 *
150 * Revision 1.10 2003/10/05 21:11:22 e_gourgoulhon
151 * - Added method std_spectral_base().
152 * - Removed method change_triad() from this file.
153 *
154 * Revision 1.9 2003/10/03 15:09:38 j_novak
155 * Improved display
156 *
157 * Revision 1.8 2003/10/03 11:21:48 j_novak
158 * More methods for the class Metric
159 *
160 * Revision 1.7 2003/10/01 11:56:31 e_gourgoulhon
161 * Corrected error: '=' replaced by '==' in two assert tests.
162 *
163 * Revision 1.6 2003/09/30 08:38:23 j_novak
164 * added a header typeinfo
165 *
166 * Revision 1.5 2003/09/29 12:52:57 j_novak
167 * Methods for changing the triad are implemented.
168 *
169 * Revision 1.4 2003/09/26 14:33:53 j_novak
170 * Arithmetic functions for the class Tensor
171 *
172 * Revision 1.3 2003/09/24 15:10:54 j_novak
173 * Suppression of the etat flag in class Tensor (still present in Scalar)
174 *
175 * Revision 1.2 2003/09/23 08:52:17 e_gourgoulhon
176 * new version
177 *
178 * Revision 1.1 2003/09/22 12:52:51 e_gourgoulhon
179 * First version: not ready yet !
180 *
181 *
182 * $Header: /cvsroot/Lorene/C++/Source/Tensor/tensor.C,v 1.43 2014/10/13 08:53:44 j_novak Exp $
183 *
184 */
185
186// Headers C
187#include <cstdlib>
188#include <cassert>
189#include <cmath>
190
191// Headers Lorene
192#include "metric.h"
193#include "utilitaires.h"
194
195 //--------------//
196 // Constructors //
197 //--------------//
198
199// Standard constructor
200// --------------------
201namespace Lorene {
202Tensor::Tensor(const Map& map, int val, const Itbl& tipe,
203 const Base_vect& triad_i)
204 : mp(&map), valence(val), triad(&triad_i), type_indice(tipe),
205 n_comp(int(pow(3., val))){
206
207 // Des verifs :
208 assert (valence >= 0) ;
209 assert (tipe.get_ndim() == 1) ;
210 assert (valence == tipe.get_dim(0)) ;
211 for (int i=0 ; i<valence ; i++)
212 assert ((tipe(i) == COV) || (tipe(i) == CON)) ;
213
214 cmp = new Scalar*[n_comp] ;
215
216 for (int i=0 ; i<n_comp ; i++)
217 cmp[i] = new Scalar(map) ;
218
219 set_der_0x0() ;
220
221}
222
223// Standard constructor with the triad passed as a pointer
224// -------------------------------------------------------
225Tensor::Tensor(const Map& map, int val, const Itbl& tipe,
226 const Base_vect* triad_i)
227 : mp(&map), valence(val), triad(triad_i), type_indice(tipe),
228 n_comp(int(pow(3., val))){
229
230 // Des verifs :
231 assert (valence >= 0) ;
232 assert (tipe.get_ndim() == 1) ;
233 assert (valence == tipe.get_dim(0)) ;
234 for (int i=0 ; i<valence ; i++)
235 assert ((tipe(i) == COV) || (tipe(i) == CON)) ;
236
237 cmp = new Scalar*[n_comp] ;
238
239 for (int i=0 ; i<n_comp ; i++)
240 cmp[i] = new Scalar(map) ;
241
242 set_der_0x0() ;
243
244}
245
246
247
248
249// Standard constructor when all the indices are of the same type
250// --------------------------------------------------------------
251Tensor::Tensor(const Map& map, int val, int tipe, const Base_vect& triad_i)
252 : mp(&map), valence(val), triad(&triad_i), type_indice(val),
253 n_comp(int(pow(3., val))){
254
255 // Des verifs :
256 assert (valence >= 0) ;
257 assert ((tipe == COV) || (tipe == CON)) ;
258
259 type_indice = tipe ;
260
261 cmp = new Scalar*[n_comp] ;
262
263 for (int i=0 ; i<n_comp ; i++)
264 cmp[i] = new Scalar(map) ;
265
266 set_der_0x0() ;
267}
268
269// Copy constructor
270// ----------------
272 mp(source.mp), valence(source.valence), triad(source.triad),
273 type_indice(source.type_indice) {
274
275 n_comp = int(pow(3., valence)) ;
276
277 cmp = new Scalar*[n_comp] ;
278 for (int i=0 ; i<n_comp ; i++) {
279
280 // the following writing takes into account the case where
281 // source belongs to a derived class of Tensor with a different
282 // storage of components :
283
285 cmp[i] = new Scalar(*source.cmp[place_source]) ;
286 }
287
288 set_der_0x0() ;
289
290}
291
292
293// Constructor from a file
294// -----------------------
296 : mp(&mapping), triad(&triad_i), type_indice(fd){
297
298 fread_be(&valence, sizeof(int), 1, fd) ;
299
300 if (valence != 0) {
302 assert( *triad_fich == *triad) ;
303 delete triad_fich ;
304 }
305 else{
306 triad = 0x0 ;
307 }
308
309 fread_be(&n_comp, sizeof(int), 1, fd) ;
310
311 cmp = new Scalar*[n_comp] ;
312 for (int i=0 ; i<n_comp ; i++)
313 cmp[i] = new Scalar(*mp, *(mp->get_mg()), fd) ;
314
315 set_der_0x0() ;
316}
317
318
319// Constructor for a scalar field: to be used by the derived
320// class {\tt Scalar}
321//-----------------------------------------------------------
322Tensor::Tensor(const Map& map) : mp(&map), valence(0), triad(0x0),
323 type_indice(0), n_comp(1) {
324
325 cmp = new Scalar*[n_comp] ;
326 cmp[0] = 0x0 ;
327
328 set_der_0x0() ;
329}
330
331
332// Constructor used by the derived classes
333// ---------------------------------------
334Tensor::Tensor (const Map& map, int val, const Itbl& tipe, int compo,
335 const Base_vect& triad_i) :
336 mp(&map), valence(val), triad(&triad_i), type_indice(tipe), n_comp(compo)
337{
338
339 // Des verifs :
340 assert (valence >= 0) ;
341 assert (tipe.get_ndim() == 1) ;
342 assert (n_comp > 0) ;
343 assert (valence == tipe.get_dim(0)) ;
344 for (int i=0 ; i<valence ; i++)
345 assert ((tipe(i) == COV) || (tipe(i) == CON)) ;
346
347 cmp = new Scalar*[n_comp] ;
348
349 for (int i=0 ; i<n_comp ; i++)
350 cmp[i] = new Scalar(map) ;
351
352 set_der_0x0() ;
353
354}
355
356// Constructor used by the derived classes when all the indices are of
357// the same type.
358// -------------------------------------------------------------------
359Tensor::Tensor (const Map& map, int val, int tipe, int compo,
360 const Base_vect& triad_i) :
361 mp(&map), valence(val), triad(&triad_i), type_indice(val), n_comp(compo)
362{
363
364 // Des verifs :
365 assert (valence >= 0) ;
366 assert (n_comp >= 0) ;
367 assert ((tipe == COV) || (tipe == CON)) ;
368
369 type_indice = tipe ;
370
371 cmp = new Scalar*[n_comp] ;
372
373 for (int i=0 ; i<n_comp ; i++)
374 cmp[i] = new Scalar(map) ;
375
376 set_der_0x0() ;
377
378}
379
380 //--------------//
381 // Destructor //
382 //--------------//
383
384
386
388
389 for (int i=0 ; i<n_comp ; i++) {
390 if (cmp[i] != 0x0)
391 delete cmp[i] ;
392 }
393 delete [] cmp ;
394}
395
396
397
398void Tensor::del_deriv() const {
399
400 for (int i=0; i<N_MET_MAX; i++)
402
403 set_der_0x0() ;
404
405}
406
408
409 for (int i=0; i<N_MET_MAX; i++)
411
412}
413
414void Tensor::del_derive_met(int j) const {
415
416 assert( (j>=0) && (j<N_MET_MAX) ) ;
417
418 if (met_depend[j] != 0x0) {
419 for (int i=0 ; i<N_TENSOR_DEPEND ; i++)
420 if (met_depend[j]->tensor_depend[i] == this)
421 met_depend[j]->tensor_depend[i] = 0x0 ;
422 if (p_derive_cov[j] != 0x0)
423 delete p_derive_cov[j] ;
424 if (p_derive_con[j] != 0x0)
425 delete p_derive_con[j] ;
426 if (p_divergence[j] != 0x0)
427 delete p_divergence[j] ;
428
430 }
431}
432
433void Tensor::set_der_met_0x0(int i) const {
434
435 assert( (i>=0) && (i<N_MET_MAX) ) ;
436 met_depend[i] = 0x0 ;
437 p_derive_cov[i] = 0x0 ;
438 p_derive_con[i] = 0x0 ;
439 p_divergence[i] = 0x0 ;
440
441}
442
444 int resu = -1 ;
445 for (int i=0; i<N_MET_MAX; i++)
446 if (met_depend[i] == &metre) {
447 resu = i ;
448 break ;
449 }
450 return resu ;
451}
452
453void Tensor::set_dependance (const Metric& met) const {
454
455 int nmet = 0 ;
456 bool deja = false ;
457 for (int i=0; i<N_MET_MAX; i++) {
458 if (met_depend[i] == &met) deja = true ;
459 if ((!deja) && (met_depend[i] != 0x0)) nmet++ ;
460 }
461 if (nmet == N_MET_MAX) {
462 cout << "Too many metrics in Tensor::set_dependances" << endl ;
463 abort() ;
464 }
465 if (!deja) {
466 int conte = 0 ;
467 while ((conte < N_TENSOR_DEPEND) && (met.tensor_depend[conte] != 0x0))
468 conte ++ ;
469
470 if (conte == N_TENSOR_DEPEND) {
471 cout << "Too many dependancies in Tensor::set_dependances " << endl ;
472 abort() ;
473 }
474 else {
475 met.tensor_depend[conte] = this ;
476 met_depend[nmet] = &met ;
477 }
478 }
479}
480
482
483 del_deriv() ;
484 for (int i=0 ; i<n_comp ; i++) {
485 cmp[i]->set_etat_qcq() ;
486 }
487}
488
490
491 del_deriv() ;
492 for (int i=0 ; i<n_comp ; i++) {
493 cmp[i]->set_etat_nondef() ;
494 }
495}
496
498
499 del_deriv() ;
500 for (int i=0 ; i<n_comp ; i++) {
501 cmp[i]->set_etat_zero() ;
502 }
503}
504
505
506// Allocates everything
507// --------------------
509
510 del_deriv() ;
511 for (int i=0 ; i<n_comp ; i++) {
512 cmp[i]->allocate_all() ;
513 }
514
515}
516
517
518
520
521 triad = &bi ;
522
523}
524
525int Tensor::position (const Itbl& idx) const {
526
527 assert (idx.get_ndim() == 1) ;
528 assert (idx.get_dim(0) == valence) ;
529
530 for (int i=0 ; i<valence ; i++)
531 assert ((idx(i)>=1) && (idx(i)<=3)) ;
532 int res = 0 ;
533 for (int i=0 ; i<valence ; i++)
534 res = 3*res+(idx(i)-1) ;
535
536 return res;
537}
538
540
541 assert ((place >= 0) && (place < n_comp)) ;
542
543 Itbl res(valence) ;
544
545 for (int i=valence-1 ; i>=0 ; i--) {
546 res.set(i) = div(place, 3).rem ;
547 place = int((place-res(i))/3) ;
548 res.set(i)++ ;
549 }
550 return res ;
551}
552
553void Tensor::operator=(const Tensor& t) {
554
555 assert (valence == t.valence) ;
556
557 triad = t.triad ;
558
559 for (int i=0 ; i<valence ; i++)
561
562 for (int i=0 ; i<n_comp ; i++) {
563 int place_t = t.position(indices(i)) ;
564 *cmp[i] = *t.cmp[place_t] ;
565 }
566
567 del_deriv() ;
568
569}
570
572
573 assert (valence == t.valence) ;
574 assert (triad == t.triad) ;
575 for (int i=0 ; i<valence ; i++)
577
578 for (int i=0 ; i<n_comp ; i++) {
579 int place_t = t.position(indices(i)) ;
580 *cmp[i] += *t.cmp[place_t] ;
581 }
582
583 del_deriv() ;
584
585}
586
588
589 assert (valence == t.valence) ;
590 assert (triad == t.triad) ;
591 for (int i=0 ; i<valence ; i++)
593
594 for (int i=0 ; i<n_comp ; i++) {
595 int place_t = t.position(indices(i)) ;
596 *cmp[i] -= *t.cmp[place_t] ;
597 }
598
599 del_deriv() ;
600
601}
602
603
604
605// Affectation d'un tenseur d'ordre 2 :
607
608 assert (valence == 2) ;
609
610 Itbl ind (valence) ;
611 ind.set(0) = ind1 ;
612 ind.set(1) = ind2 ;
613
614 int place = position(ind) ;
615
616 del_deriv() ;
617 return *cmp[place] ;
618}
619
620// Affectation d'un tenseur d'ordre 3 :
621Scalar& Tensor::set(int ind1, int ind2, int ind3) {
622
623 assert (valence == 3) ;
624
625 Itbl idx(valence) ;
626 idx.set(0) = ind1 ;
627 idx.set(1) = ind2 ;
628 idx.set(2) = ind3 ;
629 int place = position(idx) ;
630 del_deriv() ;
631
632 return *cmp[place] ;
633}
634
635
636// Affectation d'un tenseur d'ordre 4 :
637Scalar& Tensor::set(int ind1, int ind2, int ind3, int ind4) {
638
639 assert (valence == 4) ;
640
641 Itbl idx(valence) ;
642 idx.set(0) = ind1 ;
643 idx.set(1) = ind2 ;
644 idx.set(2) = ind3 ;
645 idx.set(3) = ind4 ;
646 int place = position(idx) ;
647 del_deriv() ;
648
649 return *cmp[place] ;
650}
651
652
653// Affectation cas general
655
656 assert (idx.get_ndim() == 1) ;
657 assert (idx.get_dim(0) == valence) ;
658
659 int place = position(idx) ;
660
661 del_deriv() ;
662 return *cmp[place] ;
663}
664
665// Annulation dans des domaines
667
668 annule(l, l) ;
669}
670
671void Tensor::annule(int l_min, int l_max) {
672
673 // Cas particulier: annulation globale :
674 if ( (l_min == 0) && (l_max == mp->get_mg()->get_nzone()-1) ) {
675 set_etat_zero() ;
676 return ;
677 }
678
679 // Annulation des composantes:
680 for (int i=0 ; i<n_comp ; i++) {
681 cmp[i]->annule(l_min, l_max) ;
682 }
683
684 // The derived members are no longer up to date:
685 del_deriv() ;
686
687}
688
689
691
692 // Not applicable in the nucleus nor the CED:
693 assert( mp->get_mg()->get_type_r(lrac) == FIN ) ;
694
695 int nz = mp->get_mg()->get_nzone() ;
696#ifndef NDEBUG
697 if ((2*deg+1) >= mp->get_mg()->get_nr(lrac))
698 cout << "Tensor::annule_extern_cn : \n"
699 << "WARNING!! \n"
700 << "The number of coefficients in r is too low \n"
701 << "to do a clean matching..." << endl ;
702#endif
703 // Boundary of domain lrac
704 double r_min = mp->val_r(lrac, -1., 0., 0.) ;
705 double r_max = mp->val_r(lrac, 1., 0., 0.) ;
706
707 //Definition of binomial coefficients array
708 Itbl binom(deg+1, deg+1) ;
709 binom.annule_hard() ;
710 binom.set(0,0) = 1 ;
711 for (int n=1; n<=deg; n++) {
712 binom.set(n,0) = 1 ;
713 for (int k=1; k<=n; k++)
714 binom.set(n,k) = binom(n-1, k) + binom(n-1, k-1) ;
715 }
716
717 // Coefficient of the second polynomial factor
718 Tbl coef(deg+1) ;
719 coef.set_etat_qcq() ;
720 coef.set(deg) = 1 ;
721 int sg = -1 ;
722 for (int i=deg-1; i>=0; i--) {
723
724 coef.set(i) = double(r_max*(i+1)*coef(i+1)
725 + sg*binom(deg,i)*(2*deg+1)*pow(r_min,deg-i))
726 / double(deg+i+1) ;
727 sg *= -1 ;
728 }
729
730 // Normalization to have 1 at r_min
731 double aa = coef(deg) ;
732 for (int i = deg-1; i>=0; i--)
733 aa = r_min*aa + coef(i) ;
734 aa *= pow(r_min - r_max, deg+1) ;
735 aa = 1/aa ;
736
737 Mtbl mr = mp->r ;
738 Tbl rr = mr(lrac) ;
739
740 Tbl poly(rr) ;
741 poly = coef(deg) ;
742 for (int i=deg-1; i>=0; i--)
743 poly = rr*poly + coef(i) ;
744 poly *= aa*pow((rr-r_max), deg+1) ;
745
746 Scalar rac(*mp) ;
747 rac.allocate_all() ;
748 for (int l=0; l<lrac; l++) rac.set_domain(l) = 1 ;
749 rac.set_domain(lrac) = poly ;
750 rac.annule(lrac+1,nz-1) ;
751 rac.std_spectral_base() ;
752
753 for (int ic=0; ic<n_comp; ic++) *(cmp[ic]) *= rac ;
754
755 del_deriv() ;
756}
757
758
759
760const Scalar& Tensor::operator()(int indice1, int indice2) const {
761
762 assert(valence == 2) ;
763
764 Itbl idx(2) ;
765 idx.set(0) = indice1 ;
766 idx.set(1) = indice2 ;
767 return *cmp[position(idx)] ;
768
769}
770
771const Scalar& Tensor::operator()(int indice1, int indice2, int indice3) const {
772
773 assert(valence == 3) ;
774
775 Itbl idx(3) ;
776 idx.set(0) = indice1 ;
777 idx.set(1) = indice2 ;
778 idx.set(2) = indice3 ;
779 return *cmp[position(idx)] ;
780}
781
782
784 int indice4) const {
785
786 assert(valence == 4) ;
787
788 Itbl idx(4) ;
789 idx.set(0) = indice1 ;
790 idx.set(1) = indice2 ;
791 idx.set(2) = indice3 ;
792 idx.set(3) = indice4 ;
793 return *cmp[position(idx)] ;
794}
795
796
797
798const Scalar& Tensor::operator()(const Itbl& ind) const {
799
800 assert (ind.get_ndim() == 1) ;
801 assert (ind.get_dim(0) == valence) ;
802 return *cmp[position(ind)] ;
803
804}
805
806
807// Gestion de la CED :
809
810 for (int i=0 ; i<n_comp ; i++)
812
813 del_deriv() ;
814}
815
817
818 for (int i=0 ; i<n_comp ; i++)
819 cmp[i]->inc_dzpuis(inc) ;
820
821 del_deriv() ;
822}
823
824
825// Le cout :
826ostream& operator<<(ostream& flux, const Tensor &source ) {
827
828 flux << '\n' ;
829 flux << "Lorene class : " << typeid(source).name()
830 << " Valence : " << source.valence << '\n' ;
831
832 if (source.get_triad() != 0x0) {
833 flux << "Vectorial basis (triad) on which the components are defined :"
834 << '\n' ;
835 flux << *(source.get_triad()) << '\n' ;
836 }
837
838 if (source.valence != 0) {
839 flux << "Type of the indices : " ;
840 for (int i=0 ; i<source.valence ; i++) {
841 flux << "index " << i << " : " ;
842 if (source.type_indice(i) == CON)
843 flux << " contravariant." << '\n' ;
844 else
845 flux << " covariant." << '\n' ;
846 if ( i < source.valence-1 ) flux << " " ;
847 }
848 flux << '\n' ;
849 }
850
851 for (int i=0 ; i<source.n_comp ; i++) {
852
853 if (source.valence == 0) {
854 flux <<
855 "===================== Scalar field ========================= \n" ;
856 }
857 else {
858 flux << "================ Component " ;
859 Itbl num_indices = source.indices(i) ;
860 for (int j=0 ; j<source.valence ; j++) {
861 flux << " " << num_indices(j) ;
862 }
863 flux << " ================ \n" ;
864 }
865 flux << '\n' ;
866
867 flux << *source.cmp[i] << '\n' ;
868 }
869
870 return flux ;
871}
872
873
875 double thres, int precis, ostream& ost) const {
876
877 if (comment != 0x0) {
878 ost << comment << " : " << endl ;
879 }
880
881 ost << "Lorene class : " << typeid(*this).name()
882 << " Valence : " << valence << '\n' ;
883
884 for (int ic=0; ic<n_comp; ic++) {
885
886 if (valence == 0) {
887 ost <<
888 "===================== Scalar field ========================= \n" ;
889 }
890 else {
891 ost << "================ Component " ;
893 for (int j=0 ; j<valence ; j++) {
894 ost << " " << num_indices(j) ;
895 }
896 ost << " ================ \n" ;
897 }
898 ost << '\n' ;
899
900 cmp[ic]->spectral_display(0x0, thres, precis, ost) ;
901 ost << '\n' ;
902 }
903}
904
905
906void Tensor::sauve(FILE* fd) const {
907
908 type_indice.sauve(fd) ; // type des composantes
909 fwrite_be(&valence, sizeof(int), 1, fd) ; // la valence
910
911 if (valence != 0) {
912 triad->sauve(fd) ; // Vectorial basis
913 }
914
915 fwrite_be(&n_comp, sizeof(int), 1, fd) ; // nbre composantes
916 for (int i=0 ; i<n_comp ; i++)
917 cmp[i]->sauve(fd) ;
918
919}
920
921
922
923
924
925// Sets the standard spectal bases of decomposition for each component
927
928 switch (valence) {
929
930 case 0 : {
931 cmp[0]->std_spectral_base() ;
932 break ;
933 }
934
935 case 1 : {
936 cout <<
937 "Tensor::std_spectral_base: should not be called on a Tensor"
938 << " of valence 1 but on a Vector !" << endl ;
939 abort() ;
940 break ;
941 }
942
943 case 2 : {
944
945 Base_val** bases = 0x0 ;
946 if( triad->identify() == (mp->get_bvect_cart()).identify() ) {
947 bases = mp->get_mg()->std_base_vect_cart() ;
948 }
949 else {
950 assert( triad->identify() == (mp->get_bvect_spher()).identify()) ;
951 bases = mp->get_mg()->std_base_vect_spher() ;
952 }
953
954 Itbl ind(2) ;
955 for (int i=0 ; i<n_comp ; i++) {
956 ind = indices(i) ;
957 cmp[i]->set_spectral_base( (*bases[ind(0)-1]) *
958 (*bases[ind(1)-1]) ) ;
959 }
960
961 for (int i=0 ; i<3 ; i++) {
962 delete bases[i] ;
963 }
964 delete [] bases ;
965 break ;
966
967 }
968
969
970 default : {
971
972 cout << "Tensor::std_spectral_base: the case valence = " << valence
973 << " is not treated yet !" << endl ;
974 abort() ;
975 break ;
976 }
977 }
978}
979
980// Sets the standard spectal bases of decomposition for each component (odd in the nucleus)
981
983
984 switch (valence) {
985
986 case 0 : {
988 break ;
989 }
990
991 default : {
992
993 cout << "Tensor::std_spectral_base_odd: the case valence = " << valence
994 << " is not treated yet !" << endl ;
995 abort() ;
996 break ;
997 }
998 }
999}
1000
1001
1003
1005 int j = get_place_met(metre) ;
1006 assert ((j>=0) && (j<N_MET_MAX)) ;
1007 if (p_derive_cov[j] == 0x0) {
1008 p_derive_cov[j] = metre.connect().p_derive_cov(*this) ;
1009 }
1010 return *p_derive_cov[j] ;
1011}
1012
1013
1015
1017 int j = get_place_met(metre) ;
1018 assert ((j>=0) && (j<N_MET_MAX)) ;
1019 if (p_derive_con[j] == 0x0) {
1020
1021 if (valence == 0) {
1022 p_derive_con[j] =
1023 new Vector( contract(derive_cov(metre), 0, metre.con(), 0) ) ;
1024 }
1025 else {
1026 const Tensor_sym* tsym = dynamic_cast<const Tensor_sym*>(this) ;
1027
1028 if (tsym != 0x0) { // symmetric case, preserved by derive_con
1029 const Tensor& dercov = derive_cov(metre) ;
1030 Itbl type_ind = dercov.get_index_type() ;
1031 type_ind.set(valence) = CON ;
1033 tsym->sym_index1(), tsym->sym_index2()) ;
1034
1035 *(p_derive_con[j]) = contract(dercov, valence, metre.con(), 0) ;
1036 // valence is the number of the last index of derive_cov
1037 // (the "derivation" index)
1038 }
1039 else { // general case, no symmetry
1040
1042 metre.con(), 0) ) ;
1043 // valence is the number of the last index of derive_cov
1044 // (the "derivation" index)
1045 }
1046
1047 }
1048
1049 }
1050
1051 return *p_derive_con[j] ;
1052
1053}
1054
1056
1058 int j = get_place_met(metre) ;
1059 assert ((j>=0) && (j<N_MET_MAX)) ;
1060 if (p_divergence[j] == 0x0) {
1061 p_divergence[j] = metre.connect().p_divergence(*this) ;
1062 }
1063 return *p_divergence[j] ;
1064}
1065
1067 double alpha) {
1068 if( triad->identify() == (mp->get_bvect_cart()).identify() )
1069 for (int i=0; i<n_comp; i++)
1070 cmp[i]->exponential_filter_r(lzmin, lzmax, p, alpha) ;
1071 else {
1072 cout << "Tensor::exponential_filter_r : " << endl ;
1073 cout << "Only Cartesian triad is implemented!" << endl ;
1074 cout << "Exiting..." << endl ;
1075 abort() ;
1076 }
1077}
1078
1080 double alpha) {
1081 if( triad->identify() == (mp->get_bvect_cart()).identify() )
1082 for (int i=0; i<n_comp; i++)
1083 cmp[i]->exponential_filter_ylm(lzmin, lzmax, p, alpha) ;
1084 else {
1085 cout << "Tensor::exponential_filter_ylm : " << endl ;
1086 cout << "Only Cartesian triad is implemented!" << endl ;
1087 cout << "Exiting..." << endl ;
1088 abort() ;
1089 }
1090}
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105}
Bases of the spectral expansions.
Definition base_val.h:322
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition base_vect.h:105
static Base_vect * bvect_from_file(FILE *)
Construction of a vectorial basis from a file (see sauve(FILE* ) ).
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
int position(int j) const
Gives the position in the arrays step, the_time and val corresponding to the time step j.
Definition evolution.C:273
Basic integer array class.
Definition itbl.h:122
void sauve(FILE *) const
Save in a file.
Definition itbl.C:226
Base class for coordinate mappings.
Definition map.h:670
Metric for tensor calculation.
Definition metric.h:90
Multi-domain array.
Definition mtbl.h:118
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:784
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition scalar.C:391
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition scalar.C:365
virtual void set_etat_nondef()
Sets the logical state to ETATNONDEF (undefined).
Definition scalar.C:344
virtual void std_spectral_base_odd()
Sets the spectral bases of the Valeur va to the standard odd ones for a scalar field.
Definition scalar.C:791
virtual void spectral_display(const char *comment=0x0, double threshold=1.e-7, int precision=4, ostream &ostr=cout) const
Displays the spectral coefficients and the associated basis functions.
Definition scalar.C:741
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
Definition scalar.C:797
Basic array class.
Definition tbl.h:161
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
Symmetric tensors (with respect to two of their arguments).
Definition tensor.h:1037
Tensor handling.
Definition tensor.h:288
Tensor field of valence 1.
Definition vector.h:188
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition fread_be.C:69
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition fwrite_be.C:70
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 tensor.C:1079
virtual void operator=(const Tensor &)
Assignment to a Tensor.
Definition tensor.C:553
virtual void sauve(FILE *) const
Save in a binary file.
Definition tensor.C:906
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
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition tensor.C:519
virtual void allocate_all()
Performs the memory allocation of all the elements, down to the double arrays of the Tbl s.
Definition tensor.C:508
virtual void std_spectral_base_odd()
Sets the standard odd spectal bases of decomposition for each component.
Definition tensor.C:982
Tensor * p_divergence[N_MET_MAX]
Array of pointers on the divergence of this with respect to various metrics.
Definition tensor.h:351
void annule_extern_cn(int l_0, int deg)
Performs a smooth (C^n) transition in a given domain to zero.
Definition tensor.C:690
int get_place_met(const Metric &) const
Returns the position of the pointer on metre in the array met_depend .
Definition tensor.C:443
void set_der_0x0() const
Sets the pointers on derived quantities to 0x0.
Definition tensor.C:407
virtual void spectral_display(const char *comment=0x0, double threshold=1.e-7, int precision=4, ostream &ostr=cout) const
Displays the spectral coefficients and the associated basis functions of each component.
Definition tensor.C:874
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
Tensor * p_derive_cov[N_MET_MAX]
Array of pointers on the covariant derivatives of this with respect to various metrics.
Definition tensor.h:335
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
virtual Itbl indices(int pos) const
Returns the indices of a component given by its position in the array cmp .
Definition tensor.C:539
virtual int position(const Itbl &ind) const
Returns the position in the array cmp of a component given by its indices.
Definition tensor.C:525
Tensor(const Map &map, int val, const Itbl &tipe, const Base_vect &triad_i)
Standard constructor.
Definition tensor.C:202
int valence
Valence of the tensor (0 = scalar, 1 = vector, etc...)
Definition tensor.h:298
virtual void set_etat_nondef()
Sets the logical state of all components to ETATNONDEF (undefined state).
Definition tensor.C:489
virtual void annule(int l_min, int l_max)
Sets the Tensor to zero in several domains.
Definition tensor.C:671
virtual ~Tensor()
Destructor.
Definition tensor.C:385
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:816
const Scalar & operator()(const Itbl &ind) const
Returns the value of a component (read-only version).
Definition tensor.C:798
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
int n_comp
Number of stored components, depending on the symmetry.
Definition tensor.h:312
void operator-=(const Tensor &)
-= Tensor
Definition tensor.C:587
void set_der_met_0x0(int) const
Sets all the i-th components of met_depend , p_derive_cov , etc... to 0x0.
Definition tensor.C:433
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
virtual void del_deriv() const
Deletes the derived quantities.
Definition tensor.C:398
void operator+=(const Tensor &)
+= Tensor
Definition tensor.C:571
const Tensor & divergence(const Metric &gam) const
Computes the divergence of this with respect to some metric .
Definition tensor.C:1055
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition tensor.C:481
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
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 tensor.C:1066
Tensor * p_derive_con[N_MET_MAX]
Array of pointers on the contravariant derivatives of this with respect to various metrics.
Definition tensor.h:343
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition tensor.C:497
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tensor.h:303
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition tensor.C:926
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64