KADATH
metric_harmonic.cpp
1 /*
2  Copyright 2017 Philippe Grandclement
3 
4  This file is part of Kadath.
5 
6  Kadath is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  Kadath is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with Kadath. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #include "space.hpp"
21 #include "tensor.hpp"
22 #include "metric.hpp"
23 #include "term_eq.hpp"
24 #include "scalar.hpp"
25 #include "tensor_impl.hpp"
26 #include "system_of_eqs.hpp"
27 #include "metric_tensor.hpp"
28 namespace Kadath {
29 Metric_harmonic::Metric_harmonic (Metric_tensor& met) :
30  Metric_general(met) {
31  type_tensor = met.get_type() ;
32 }
33 
34 Metric_harmonic::Metric_harmonic (const Metric_harmonic& so) :
35  Metric_general (so) {
36 }
37 
38 Metric_harmonic::~Metric_harmonic() {
39 }
40 
41 void Metric_harmonic::compute_det_cov (int dd) const {
42 
43  // Need that
44  if (p_met_cov[dd]==0x0)
45  compute_cov(dd) ;
46 
47  bool doder = (p_met_cov[dd]->der_t==0x0) ? false : true ;
48 
49  Scalar res_val (espace) ;
50  Scalar res_der (espace) ;
51 
52  Val_domain cmpval (espace.get_domain(dd)) ;
53  cmpval = 0 ;
54 
55  for (int i=1 ; i<=espace.get_ndim() ; i++)
56  for (int j=1 ; j<=espace.get_ndim() ; j++)
57  cmpval += (*p_met_cov[dd]->val_t)(1,1)(dd)*(*p_met_cov[dd]->val_t)(2,2)(dd)*(*p_met_cov[dd]->val_t)(3,3)(dd)
58  + (*p_met_cov[dd]->val_t)(1,2)(dd)*(*p_met_cov[dd]->val_t)(2,3)(dd)*(*p_met_cov[dd]->val_t)(1,3)(dd)
59  + (*p_met_cov[dd]->val_t)(1,3)(dd)*(*p_met_cov[dd]->val_t)(1,2)(dd)*(*p_met_cov[dd]->val_t)(2,3)(dd)
60  - (*p_met_cov[dd]->val_t)(1,3)(dd)*(*p_met_cov[dd]->val_t)(1,3)(dd)*(*p_met_cov[dd]->val_t)(2,2)(dd)
61  - (*p_met_cov[dd]->val_t)(2,3)(dd)*(*p_met_cov[dd]->val_t)(2,3)(dd)*(*p_met_cov[dd]->val_t)(1,1)(dd)
62  - (*p_met_cov[dd]->val_t)(1,2)(dd)*(*p_met_cov[dd]->val_t)(1,2)(dd)*(*p_met_cov[dd]->val_t)(3,3)(dd) ;
63  res_val.set_domain(dd) = cmpval ;
64 
65  if (doder) {
66  Val_domain cmpder (espace.get_domain(dd)) ;
67  cmpder = 0 ;
68 
69  for (int i=1 ; i<=espace.get_ndim() ; i++)
70  for (int j=1 ; j<=espace.get_ndim() ; j++)
71  cmpder += (*p_met_cov[dd]->der_t)(1,1)(dd)*(*p_met_cov[dd]->val_t)(2,2)(dd)*(*p_met_cov[dd]->val_t)(3,3)(dd)
72  + (*p_met_cov[dd]->der_t)(1,2)(dd)*(*p_met_cov[dd]->val_t)(2,3)(dd)*(*p_met_cov[dd]->val_t)(1,3)(dd)
73  + (*p_met_cov[dd]->der_t)(1,3)(dd)*(*p_met_cov[dd]->val_t)(1,2)(dd)*(*p_met_cov[dd]->val_t)(2,3)(dd)
74  - (*p_met_cov[dd]->der_t)(1,3)(dd)*(*p_met_cov[dd]->val_t)(1,3)(dd)*(*p_met_cov[dd]->val_t)(2,2)(dd)
75  - (*p_met_cov[dd]->der_t)(2,3)(dd)*(*p_met_cov[dd]->val_t)(2,3)(dd)*(*p_met_cov[dd]->val_t)(1,1)(dd)
76  - (*p_met_cov[dd]->der_t)(1,2)(dd)*(*p_met_cov[dd]->val_t)(1,2)(dd)*(*p_met_cov[dd]->val_t)(3,3)(dd)
77  + (*p_met_cov[dd]->val_t)(1,1)(dd)*(*p_met_cov[dd]->der_t)(2,2)(dd)*(*p_met_cov[dd]->val_t)(3,3)(dd)
78  + (*p_met_cov[dd]->val_t)(1,2)(dd)*(*p_met_cov[dd]->der_t)(2,3)(dd)*(*p_met_cov[dd]->val_t)(1,3)(dd)
79  + (*p_met_cov[dd]->val_t)(1,3)(dd)*(*p_met_cov[dd]->der_t)(1,2)(dd)*(*p_met_cov[dd]->val_t)(2,3)(dd)
80  - (*p_met_cov[dd]->val_t)(1,3)(dd)*(*p_met_cov[dd]->der_t)(1,3)(dd)*(*p_met_cov[dd]->val_t)(2,2)(dd)
81  - (*p_met_cov[dd]->val_t)(2,3)(dd)*(*p_met_cov[dd]->der_t)(2,3)(dd)*(*p_met_cov[dd]->val_t)(1,1)(dd)
82  - (*p_met_cov[dd]->val_t)(1,2)(dd)*(*p_met_cov[dd]->der_t)(1,2)(dd)*(*p_met_cov[dd]->val_t)(3,3)(dd)
83  + (*p_met_cov[dd]->val_t)(1,1)(dd)*(*p_met_cov[dd]->val_t)(2,2)(dd)*(*p_met_cov[dd]->der_t)(3,3)(dd)
84  + (*p_met_cov[dd]->val_t)(1,2)(dd)*(*p_met_cov[dd]->val_t)(2,3)(dd)*(*p_met_cov[dd]->der_t)(1,3)(dd)
85  + (*p_met_cov[dd]->val_t)(1,3)(dd)*(*p_met_cov[dd]->val_t)(1,2)(dd)*(*p_met_cov[dd]->der_t)(2,3)(dd)
86  - (*p_met_cov[dd]->val_t)(1,3)(dd)*(*p_met_cov[dd]->val_t)(1,3)(dd)*(*p_met_cov[dd]->der_t)(2,2)(dd)
87  - (*p_met_cov[dd]->val_t)(2,3)(dd)*(*p_met_cov[dd]->val_t)(2,3)(dd)*(*p_met_cov[dd]->der_t)(1,1)(dd)
88  - (*p_met_cov[dd]->val_t)(1,2)(dd)*(*p_met_cov[dd]->val_t)(1,2)(dd)*(*p_met_cov[dd]->der_t)(3,3)(dd) ;
89  res_der.set_domain(dd) = cmpder ;
90  }
91 
92  if (!doder) {
93  if (p_det_cov[dd]==0x0)
94  p_det_cov[dd] = new Term_eq(dd, res_val) ;
95  else
96  *p_det_cov[dd] = Term_eq (dd, res_val) ;
97  }
98  else {
99  if (p_det_cov[dd]==0x0)
100  p_det_cov[dd] = new Term_eq(dd, res_val, res_der) ;
101  else
102  *p_det_cov[dd] = Term_eq(dd, res_val, res_der) ;
103  }
104 }
105 
106 
107 void Metric_harmonic::compute_ricci_tensor (int dd) const {
108 
109  // Need christoffels
110  if (p_det_cov[dd]==0x0)
111  compute_det_cov(dd) ;
112 
113  if (p_christo[dd]==0x0)
114  compute_christo(dd) ;
115  if (p_dirac[dd]==0x0)
116  compute_dirac(dd) ;
117 
118  // Need christoffels
119  if (p_christo[dd]==0x0)
120  compute_christo(dd) ;
121 
122  Tensor res_val (espace, 2, COV, basis) ;
123  Tensor res_der (espace, 2, COV, basis) ;
124 
125  Term_eq der_cov (fmet.derive (COV, ' ', (*p_met_cov[dd]))) ;
126  Term_eq der_con (fmet.derive (COV, ' ', (*p_met_con[dd]))) ;
127  Term_eq dder_cov (fmet.derive (COV, 'u', der_cov)) ;
128 
129  bool doder = ((der_cov.der_t==0x0) || (der_con.der_t==0x0) || (dder_cov.der_t==0x0)) ? false : true ;
130 
131  Index pos (res_val) ;
132  do {
133  Val_domain cmpval (espace.get_domain(dd)) ;
134  cmpval = 0 ;
135 
136  for (int k=1 ; k<=espace.get_ndim() ; k++)
137  for (int l=1 ; l<=espace.get_ndim() ; l++) {
138 
139  cmpval += -0.5 * (*p_met_con[dd]->val_t)(k,l)(dd)*(*dder_cov.val_t)(k,l, pos(0)+1, pos(1)+1)(dd)
140  - 0.5*(*der_con.val_t)(pos(1)+1, k,l)(dd)*(*der_cov.val_t)(k, pos(0)+1, l)(dd)
141  - 0.5*(*der_con.val_t)(pos(0)+1, k,l)(dd)*(*der_cov.val_t)(k, pos(1)+1, l)(dd)
142  + 0.5*(*der_con.val_t)(pos(0)+1, k,l)(dd)*(*der_cov.val_t)(pos(1)+1, k, l)(dd)
143  + (*p_christo[dd]->val_t)(pos(0)+1, pos(1)+1, k)(dd)*(*p_christo[dd]->val_t)(k, l, l)(dd)
144  - (*p_christo[dd]->val_t)(pos(0)+1, k, l )(dd)*(*p_christo[dd]->val_t)(pos(1)+1, l, k)(dd) ;
145 
146  for (int m=1 ; m<=espace.get_ndim() ; m++) {
147  cmpval += (*p_met_cov[dd]->val_t)(m,l)(dd)* (*der_con.val_t)(k, k,l)(dd) * (*p_christo[dd]->val_t)(pos(0)+1, pos(1)+1, m)(dd)
148  - (*p_met_cov[dd]->val_t)(m,l)(dd) * (*der_con.val_t)(pos(1)+1, k,l)(dd) * (*p_christo[dd]->val_t)(pos(0)+1, k, m)(dd) ;
149 
150 
151  }
152  }
153  res_val.set(pos).set_domain(dd) = cmpval ;
154 
155  if (doder) {
156  Val_domain cmpder (espace.get_domain(dd)) ;
157  cmpder = 0 ;
158 
159  for (int k=1 ; k<=espace.get_ndim() ; k++)
160  for (int l=1 ; l<=espace.get_ndim() ; l++) {
161 
162  cmpder += -0.5 * (*p_met_con[dd]->der_t)(k,l)(dd)*(*dder_cov.val_t)(k,l, pos(0)+1, pos(1)+1)(dd)
163  -0.5 * (*p_met_con[dd]->val_t)(k,l)(dd)*(*dder_cov.der_t)(k,l, pos(0)+1, pos(1)+1)(dd)
164  - 0.5*(*der_con.der_t)(pos(1)+1, k,l)(dd)*(*der_cov.val_t)(k, pos(0)+1, l)(dd)
165  - 0.5*(*der_con.val_t)(pos(1)+1, k,l)(dd)*(*der_cov.der_t)(k, pos(0)+1, l)(dd)
166  - 0.5*(*der_con.der_t)(pos(0)+1, k,l)(dd)*(*der_cov.val_t)(k, pos(1)+1, l)(dd)
167  - 0.5*(*der_con.val_t)(pos(0)+1, k,l)(dd)*(*der_cov.der_t)(k, pos(1)+1, l)(dd)
168  + 0.5*(*der_con.der_t)(pos(0)+1, k,l)(dd)*(*der_cov.val_t)(pos(1)+1, k, l)(dd)
169  + 0.5*(*der_con.val_t)(pos(0)+1, k,l)(dd)*(*der_cov.der_t)(pos(1)+1, k, l)(dd)
170  + (*p_christo[dd]->der_t)(pos(0)+1, pos(1)+1, k)(dd)*(*p_christo[dd]->val_t)(k, l, l)(dd)
171  + (*p_christo[dd]->val_t)(pos(0)+1, pos(1)+1, k)(dd)*(*p_christo[dd]->der_t)(k, l, l)(dd)
172  - (*p_christo[dd]->der_t)(pos(0)+1, k, l )(dd)*(*p_christo[dd]->val_t)(pos(1)+1, l, k)(dd) ;
173  - (*p_christo[dd]->val_t)(pos(0)+1, k, l )(dd)*(*p_christo[dd]->der_t)(pos(1)+1, l, k)(dd) ;
174 
175  for (int m=1 ; m<=espace.get_ndim() ; m++) {
176 
177  cmpder += (*p_met_cov[dd]->der_t)(m,l)(dd)* (*der_con.val_t)(k, k,l)(dd) * (*p_christo[dd]->val_t)(pos(0)+1, pos(1)+1, m)(dd)
178  + (*p_met_cov[dd]->val_t)(m,l)(dd)* (*der_con.der_t)(k, k,l)(dd) * (*p_christo[dd]->val_t)(pos(0)+1, pos(1)+1, m)(dd)
179  + (*p_met_cov[dd]->val_t)(m,l)(dd)* (*der_con.val_t)(k, k,l)(dd) * (*p_christo[dd]->der_t)(pos(0)+1, pos(1)+1, m)(dd)
180  - (*p_met_cov[dd]->der_t)(m,l)(dd) * (*der_con.val_t)(pos(1)+1, k,l)(dd) * (*p_christo[dd]->val_t)(pos(0)+1, k, m)(dd)
181  - (*p_met_cov[dd]->val_t)(m,l)(dd) * (*der_con.der_t)(pos(1)+1, k,l)(dd) * (*p_christo[dd]->val_t)(pos(0)+1, k, m)(dd)
182  - (*p_met_cov[dd]->val_t)(m,l)(dd) * (*der_con.val_t)(pos(1)+1, k,l)(dd) * (*p_christo[dd]->der_t)(pos(0)+1, k, m)(dd) ;
183  }
184  }
185 
186 
187  res_der.set(pos).set_domain(dd) = cmpder ;
188  }
189  }
190 
191  while (pos.inc()) ;
192 
193  if (!doder) {
194  if (p_ricci_tensor[dd]==0x0)
195  p_ricci_tensor[dd] = new Term_eq(dd, res_val) ;
196  else
197  *p_ricci_tensor[dd] = Term_eq (dd, res_val) ;
198  }
199  else {
200  if (p_ricci_tensor[dd]==0x0)
201  p_ricci_tensor[dd] = new Term_eq(dd, res_val, res_der) ;
202  else
203  *p_ricci_tensor[dd] = Term_eq(dd, res_val, res_der) ;
204  }
205 
206 }
207 
208 
209 void Metric_harmonic::compute_ricci_scalar (int dd) const {
210 
211  // Need that
212  if (p_met_con[dd]==0x0)
213  compute_con(dd) ;
214  if (p_ricci_tensor[dd]==0x0)
215  compute_ricci_tensor(dd) ;
216 
217  bool doder = ((p_met_con[dd]->der_t==0x0) || (p_ricci_tensor[dd]->der_t==0x0)) ? false : true ;
218  Scalar res_val (espace) ;
219  Scalar res_der (espace) ;
220 
221  Val_domain cmpval (espace.get_domain(dd)) ;
222  cmpval = 0 ;
223 
224  for (int i=1 ; i<=espace.get_ndim() ; i++)
225  for (int j=1 ; j<=espace.get_ndim() ; j++)
226  cmpval += (*p_met_con[dd]->val_t)(i,j)(dd)*(*p_ricci_tensor[dd]->val_t)(i,j)(dd) ;
227  res_val.set_domain(dd) = cmpval ;
228 
229  if (doder) {
230  Val_domain cmpder (espace.get_domain(dd)) ;
231  cmpder = 0 ;
232 
233  for (int i=1 ; i<=espace.get_ndim() ; i++)
234  for (int j=1 ; j<=espace.get_ndim() ; j++)
235  cmpder += (*p_met_con[dd]->val_t)(i,j)(dd)*(*p_ricci_tensor[dd]->der_t)(i,j)(dd)
236  + (*p_met_con[dd]->der_t)(i,j)(dd)*(*p_ricci_tensor[dd]->val_t)(i,j)(dd) ;
237  res_der.set_domain(dd) = cmpder ;
238  }
239 
240  if (!doder) {
241  if (p_ricci_scalar[dd]==0x0)
242  p_ricci_scalar[dd] = new Term_eq(dd, res_val) ;
243  else
244  *p_ricci_scalar[dd] = Term_eq (dd, res_val) ;
245  }
246  else {
247  if (p_ricci_scalar[dd]==0x0)
248  p_ricci_scalar[dd] = new Term_eq(dd, res_val, res_der) ;
249  else
250  *p_ricci_scalar[dd] = Term_eq(dd, res_val, res_der) ;
251  }
252 }}
Kadath::Index
Class that gives the position inside a multi-dimensional Array.
Definition: index.hpp:38
Kadath::Index::inc
bool inc(int increm, int var=0)
Increments the position of the Index.
Definition: index.hpp:99
Kadath::Metric_flat::derive
virtual Term_eq derive(int, char, const Term_eq &) const
Computes the covariant derivative of a Term_eq (assumes Cartesian basis of decomposition).
Definition: metric_flat.cpp:422
Kadath::Metric_general
Class to deal with arbitrary type of metric.
Definition: metric.hpp:379
Kadath::Metric_general::compute_dirac
virtual void compute_dirac(int) const
Computes the Dirac gauge term, in a given Domain.
Definition: metric_general.cpp:756
Kadath::Metric_general::fmet
Metric_flat fmet
Associated flat metric.
Definition: metric.hpp:384
Kadath::Metric_general::compute_christo
virtual void compute_christo(int) const
Computes the Christoffel symbols, in a given Domain.
Definition: metric_general.cpp:428
Kadath::Metric_general::compute_cov
virtual void compute_cov(int) const
Computes the covariant representation, in a given Domain.
Definition: metric_general.cpp:49
Kadath::Metric_general::basis
const Base_tensor & basis
The tensorial basis used.
Definition: metric.hpp:383
Kadath::Metric_general::compute_con
virtual void compute_con(int) const
Computes the contravariant representation, in a given Domain.
Definition: metric_general.cpp:238
Kadath::Metric_harmonic
Class to deal with a metric in the spatial harmonic gauge.
Definition: metric.hpp:535
Kadath::Metric_harmonic::compute_det_cov
virtual void compute_det_cov(int) const
Computes the determinant of the covariant representation, in a given Domain.
Definition: metric_harmonic.cpp:41
Kadath::Metric_harmonic::Metric_harmonic
Metric_harmonic(Metric_tensor &)
Constructor from a Metric_tensor.
Definition: metric_harmonic.cpp:29
Kadath::Metric_harmonic::compute_ricci_scalar
virtual void compute_ricci_scalar(int) const
Computes the Ricci scalar, in a given Domain.
Definition: metric_harmonic.cpp:209
Kadath::Metric_harmonic::compute_ricci_tensor
virtual void compute_ricci_tensor(int) const
Computes the Ricci tensor, in a given Domain.
Definition: metric_harmonic.cpp:107
Kadath::Metric_tensor
Particular type of Tensor, dedicated to the desription of metrics.
Definition: metric_tensor.hpp:40
Kadath::Metric_tensor::get_type
int get_type() const
Definition: metric_tensor.hpp:81
Kadath::Metric::type_tensor
int type_tensor
States if one works in the CON or COV representation.
Definition: metric.hpp:86
Kadath::Metric::p_ricci_tensor
MMPtr_array< Term_eq > p_ricci_tensor
Array of pointers on various Term_eq.
Definition: metric.hpp:70
Kadath::Metric::espace
const Space & espace
The associated Space.
Definition: metric.hpp:42
Kadath::Metric::p_met_cov
MMPtr_array< Term_eq > p_met_cov
Array of pointers on various Term_eq.
Definition: metric.hpp:50
Kadath::Metric::p_met_con
MMPtr_array< Term_eq > p_met_con
Array of pointers on various Term_eq.
Definition: metric.hpp:55
Kadath::Metric::p_dirac
MMPtr_array< Term_eq > p_dirac
Array of pointers on various Term_eq.
Definition: metric.hpp:80
Kadath::Metric::p_christo
MMPtr_array< Term_eq > p_christo
Array of pointers on various Term_eq.
Definition: metric.hpp:60
Kadath::Metric::p_ricci_scalar
MMPtr_array< Term_eq > p_ricci_scalar
Array of pointers on various Term_eq.
Definition: metric.hpp:75
Kadath::Metric::p_det_cov
MMPtr_array< Term_eq > p_det_cov
Array of pointers on various Term_eq.
Definition: metric.hpp:85
Kadath::Scalar
The class Scalar does not really implements scalars in the mathematical sense but rather tensorial co...
Definition: scalar.hpp:67
Kadath::Scalar::set_domain
Val_domain & set_domain(int)
Read/write of a particular Val_domain.
Definition: scalar.hpp:555
Kadath::Space::get_domain
const Domain * get_domain(int i) const
returns a pointer on the domain.
Definition: space.hpp:1385
Kadath::Space::get_ndim
int get_ndim() const
Returns the number of dimensions.
Definition: space.hpp:1373
Kadath::Tensor
Tensor handling.
Definition: tensor.hpp:149
Kadath::Tensor::set
Scalar & set(const Array< int > &ind)
Returns the value of a component (read/write version).
Definition: tensor_impl.hpp:91
Kadath::Term_eq
This class is intended to describe the manage objects appearing in the equations.
Definition: term_eq.hpp:62
Kadath::Term_eq::der_t
Tensor * der_t
Pointer on the variation, if the Term_eq is a Tensor.
Definition: term_eq.hpp:69
Kadath::Term_eq::val_t
Tensor * val_t
Pointer on the value, if the Term_eq is a Tensor.
Definition: term_eq.hpp:68
Kadath::Val_domain
Class for storing the basis of decompositions of a field and its values on both the configuration and...
Definition: val_domain.hpp:69