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