KADATH
domain_shell_for_metric.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 "headcpp.hpp"
21 
22 #include "utilities.hpp"
23 #include "spheric.hpp"
24 #include "val_domain.hpp"
25 #include "scalar.hpp"
26 #include "tensor_impl.hpp"
27 #include "term_eq.hpp"
28 #include "metric.hpp"
29 namespace Kadath {
30 Term_eq Domain_shell::derive_flat_spher (int type_der, char ind_der, const Term_eq& so, const Metric* manipulator) const {
31 
32  assert ((type_der==COV) || (type_der==CON)) ;
33  int val_res = so.val_t->get_valence() + 1 ;
34 
35  bool donames = (ind_der==' ') ? false : true ;
36  bool need_sum = false ;
37 
38  Array<int> type_ind (val_res) ;
39  type_ind.set(0) = COV ;
40  for (int i=1 ; i<val_res ; i++)
41  type_ind.set(i) = so.val_t->get_index_type(i-1) ;
42 
43  if (donames) {
44  if (so.val_t->get_valence()>0)
45  assert (so.val_t->is_name_affected()) ;
46  // Need for summation ?
47  for (int i=1 ; i<val_res ; i++)
48  if (ind_der== so.val_t->get_name_ind()[i-1])
49  need_sum = true ;
50  }
51 
52  // Tensor for val
53  Base_tensor basis (so.val_t->get_space(), SPHERICAL_BASIS) ;
54  Tensor auxi_val (so.val_t->get_space(), val_res, type_ind, basis) ;
55 
56  if (donames) {
57  // Set the names of the indices :
58  auxi_val.set_name_affected() ;
59  auxi_val.set_name_ind(0, ind_der) ;
60  for (int i=1 ; i<val_res ; i++)
61  auxi_val.set_name_ind(i, so.val_t->get_name_ind()[i-1]) ;
62  }
63 
64  // Part derivative
65  //Loop on the components :
66  Index pos_auxi_bis(auxi_val) ;
67  Index pos_so_bis (*so.val_t) ;
68  do {
69  for (int i=0 ; i<val_res-1 ; i++)
70  pos_so_bis.set(i) = pos_auxi_bis(i+1) ;
71  switch (pos_auxi_bis(0)) {
72  case 0 :
73  // d/dr :
74  auxi_val.set(pos_auxi_bis).set_domain(num_dom) = (*so.val_t)(pos_so_bis)(num_dom).der_r() ;
75  break ;
76  case 1 :
77  // 1/r dtheta
78  auxi_val.set(pos_auxi_bis).set_domain(num_dom) = (*so.val_t)(pos_so_bis)(num_dom).der_var(2).div_r() ;
79  break ;
80  case 2 :
81  // 1/r sint d/dphi
82  auxi_val.set(pos_auxi_bis).set_domain(num_dom) = (*so.val_t)(pos_so_bis)(num_dom).der_var(3).div_r().div_sin_theta() ;
83  break ;
84  default :
85  cerr << "Bad indice in Domain_shell::derive_flat_spher" << endl ;
86  abort() ;
87  }
88  }
89  while (pos_auxi_bis.inc()) ;
90 
91  // Loop indice summation on connection symbols
92  for (int ind_sum=0 ; ind_sum<val_res-1 ; ind_sum++) {
93 
94  //Loop on the components :
95  Index pos_auxi(auxi_val) ;
96  Index pos_so (*so.val_t) ;
97 
98  do {
99  for (int i=0 ; i<val_res-1 ; i++)
100  pos_so.set(i) = pos_auxi(i+1) ;
101  // Different cases of the derivative index :
102  switch (pos_auxi(0)) {
103  case 0 :
104  // Dr nothing
105  break ;
106  case 1 :
107  // Dtheta
108  // Different cases of the source index
109  switch (pos_auxi(ind_sum+1)) {
110  case 0 :
111  //Dtheta S_r
112  pos_so.set(ind_sum) = 1 ;
113  auxi_val.set(pos_auxi).set_domain(num_dom) -= (*so.val_t)(pos_so)(num_dom).div_r() ;
114  break ;
115  case 1 :
116  // Dtheta S_theta
117  pos_so.set(ind_sum) = 0 ;
118  auxi_val.set(pos_auxi).set_domain(num_dom) += (*so.val_t)(pos_so)(num_dom).div_r() ;
119  break ;
120  case 2 :
121  //Dtheta S_phi
122  break ;
123  default :
124  cerr << "Bad indice in Domain_shell::derive_flat_spher" << endl ;
125  abort() ;
126  }
127  break ;
128  case 2 :
129  // Dphi
130  // Different cases of the source index
131  switch (pos_auxi(ind_sum+1)) {
132  case 0 :
133  //Dphi S_r
134  pos_so.set(ind_sum) = 2 ;
135  auxi_val.set(pos_auxi).set_domain(num_dom) -= (*so.val_t)(pos_so)(num_dom).div_r() ;
136  break ;
137  case 1 :
138  // Dphi S_theta
139  pos_so.set(ind_sum) = 2 ;
140  auxi_val.set(pos_auxi).set_domain(num_dom) -= (*so.val_t)(pos_so)(num_dom).div_r().mult_cos_theta().div_sin_theta() ;
141  break ;
142  case 2 :
143  //Dphi S_phi
144  pos_so.set(ind_sum) = 0 ;
145  auxi_val.set(pos_auxi).set_domain(num_dom) += (*so.val_t)(pos_so)(num_dom).div_r() ;
146  pos_so.set(ind_sum) = 1 ;
147  auxi_val.set(pos_auxi).set_domain(num_dom) += (*so.val_t)(pos_so)(num_dom).div_r().mult_cos_theta().div_sin_theta() ;
148  break ;
149  default :
150  cerr << "Bad indice in Domain_shell::derive_flat_spher" << endl ;
151  abort() ;
152  }
153  break ;
154  default :
155  cerr << "Bad indice in Domain_shell::derive_flat_spher" << endl ;
156  abort() ;
157  }
158  }
159  while (pos_auxi.inc()) ;
160  }
161 
162  if (so.der_t==0x0) {
163  // No need for derivative :
164  Term_eq auxi (num_dom, auxi_val) ;
165  // If derive contravariant : manipulate first indice :
166  if (type_der==CON)
167  manipulator->manipulate_ind (auxi, 0) ;
168 
169  if (!need_sum)
170  return auxi ;
171  else
172  return Term_eq (num_dom, auxi.val_t->do_summation_one_dom(num_dom)) ;
173  }
174  else {
175  // Need to compute the derivative :
176  // Tensor for der
177  Tensor auxi_der (so.val_t->get_space(), val_res, type_ind, basis) ;
178 
179  if (donames) {
180  // Set the names of the indices :
181  auxi_der.set_name_affected() ;
182  auxi_der.set_name_ind(0, ind_der) ;
183  for (int i=1 ; i<val_res ; i++)
184  auxi_der.set_name_ind(i, so.der_t->get_name_ind()[i-1]) ;
185  }
186 
187  // Part derivative
188  //Loop on the components :
189  Index pos_auxi_der_bis(auxi_der) ;
190  do {
191  for (int i=0 ; i<val_res-1 ; i++)
192  pos_so_bis.set(i) = pos_auxi_der_bis(i+1) ;
193  switch (pos_auxi_der_bis(0)) {
194  case 0 :
195  // d/dr :
196  auxi_der.set(pos_auxi_der_bis).set_domain(num_dom) = (*so.der_t)(pos_so_bis)(num_dom).der_r() ;
197  break ;
198  case 1 :
199  // 1/r dtheta
200  auxi_der.set(pos_auxi_der_bis).set_domain(num_dom) = (*so.der_t)(pos_so_bis)(num_dom).der_var(2).div_r() ;
201  break ;
202  case 2 :
203  // 1/r sint d/dphi
204  auxi_der.set(pos_auxi_der_bis).set_domain(num_dom) = (*so.der_t)(pos_so_bis)(num_dom).der_var(3).div_r().div_sin_theta() ;
205  break ;
206  default :
207  cerr << "Bad indice in Domain_shell::derive_flat_spher" << endl ;
208  abort() ;
209  }
210  }
211  while (pos_auxi_der_bis.inc()) ;
212 
213  // Loop indice summation on connection symbols
214  for (int ind_sum=0 ; ind_sum<val_res-1 ; ind_sum++) {
215 
216  //Loop on the components :
217  Index pos_auxi_der(auxi_val) ;
218  Index pos_so (*so.val_t) ;
219 
220  do {
221  for (int i=0 ; i<val_res-1 ; i++)
222  pos_so.set(i) = pos_auxi_der(i+1) ;
223  // Different cases of the derivative index :
224  switch (pos_auxi_der(0)) {
225  case 0 :
226  // Dr nothing
227  break ;
228  case 1 :
229  // Dtheta
230  // Different cases of the source index
231  switch (pos_auxi_der(ind_sum+1)) {
232  case 0 :
233  //Dtheta S_r
234  pos_so.set(ind_sum) = 1 ;
235  auxi_der.set(pos_auxi_der).set_domain(num_dom) -= (*so.der_t)(pos_so)(num_dom).div_r() ;
236  break ;
237  case 1 :
238  // Dtheta S_theta
239  pos_so.set(ind_sum) = 0 ;
240  auxi_der.set(pos_auxi_der).set_domain(num_dom) += (*so.der_t)(pos_so)(num_dom).div_r() ;
241  break ;
242  case 2 :
243  //Dtheta S_phi
244  break ;
245  default :
246  cerr << "Bad indice in Domain_shell::derive_flat_spher" << endl ;
247  abort() ;
248  }
249  break ;
250  case 2 :
251  // Dphi
252  // Different cases of the source index
253  switch (pos_auxi_der(ind_sum+1)) {
254  case 0 :
255  //Dphi S_r
256  pos_so.set(ind_sum) = 2 ;
257  auxi_der.set(pos_auxi_der).set_domain(num_dom) -= (*so.der_t)(pos_so)(num_dom).div_r() ;
258  break ;
259  case 1 :
260  // Dphi S_theta
261  pos_so.set(ind_sum) = 2 ;
262  auxi_der.set(pos_auxi_der).set_domain(num_dom) -= (*so.der_t)(pos_so)(num_dom).div_r().mult_cos_theta().div_sin_theta() ;
263  break ;
264  case 2 :
265  //Dphi S_phi
266  pos_so.set(ind_sum) = 0 ;
267  auxi_der.set(pos_auxi_der).set_domain(num_dom) += (*so.der_t)(pos_so)(num_dom).div_r() ;
268  pos_so.set(ind_sum) = 1 ;
269  auxi_der.set(pos_auxi_der).set_domain(num_dom) += (*so.der_t)(pos_so)(num_dom).div_r().mult_cos_theta().div_sin_theta() ;
270  break ;
271  default :
272  cerr << "Bad indice in Domain_shell::derive_flat_spher" << endl ;
273  abort() ;
274  }
275  break ;
276  default :
277  cerr << "Bad indice in Domain_shell::derive_flat_spher" << endl ;
278  abort() ;
279  }
280  }
281  while (pos_auxi_der.inc()) ;
282  }
283 
284 
285  // Need for derivative :
286  Term_eq auxi (num_dom, auxi_val, auxi_der) ;
287  // If derive contravariant : manipulate first indice :
288  if (type_der==CON)
289  manipulator->manipulate_ind (auxi, 0) ;
290  if (!need_sum)
291  return auxi ;
292  else
293  return Term_eq (num_dom, auxi.val_t->do_summation_one_dom(num_dom),
294  auxi.der_t->do_summation_one_dom(num_dom)) ;
295  }
296 }
297 
298 
299 Term_eq Domain_shell::derive_flat_cart (int type_der, char ind_der, const Term_eq& so, const Metric* manipulator) const {
300 
301  bool doder = (so.der_t==0x0) ? false : true ;
302  if ((type_der==CON) && (manipulator->p_met_con[num_dom]->der_t==0x0))
303  doder = false ;
304 
305  assert ((type_der==COV) || (type_der==CON)) ;
306  int val_res = so.val_t->get_valence() + 1 ;
307 
308  bool doname = true ;
309  if (so.val_t->get_valence()>0)
310  if (!so.val_t->is_name_affected())
311  doname = false;
312 
313  Array<int> type_ind (val_res) ;
314  type_ind.set(0) = COV ;
315  for (int i=1 ; i<val_res ; i++)
316  type_ind.set(i) = so.val_t->get_index_type(i-1) ;
317 
318  // Need for summation ?
319  bool need_sum = false ;
320  if (doname)
321  for (int i=1 ; i<val_res ; i++)
322  if (ind_der== so.val_t->get_name_ind()[i-1])
323  need_sum = true ;
324 
325  // Tensor for val
326  Base_tensor basis (so.val_t->get_space(), CARTESIAN_BASIS) ;
327  Tensor auxi_val (so.val_t->get_space(), val_res, type_ind, basis) ;
328 
329 
330  // Set the names of the indices :
331  if (doname) {
332  auxi_val.set_name_affected() ;
333  auxi_val.set_name_ind(0, ind_der) ;
334  for (int i=1 ; i<val_res ; i++)
335  auxi_val.set_name_ind(i, so.val_t->get_name_ind()[i-1]) ;
336  }
337 
338  //Loop on the components :
339  Index pos_auxi(auxi_val) ;
340  Index pos_so (*so.val_t) ;
341  do {
342  for (int i=0 ; i<val_res-1 ; i++)
343  pos_so.set(i) = pos_auxi(i+1) ;
344  auxi_val.set(pos_auxi).set_domain(num_dom) = (*so.val_t)(pos_so)(num_dom).der_abs(pos_auxi(0)+1) ;
345  }
346  while (pos_auxi.inc()) ;
347 
348  if (!doder) {
349  // No need for derivative :
350  Term_eq auxi (num_dom, auxi_val) ;
351  // If derive contravariant : manipulate first indice :
352  if (type_der==CON)
353  manipulator->manipulate_ind (auxi, 0) ;
354 
355  if (!need_sum)
356  return auxi ;
357  else
358  return Term_eq (num_dom, auxi.val_t->do_summation_one_dom(num_dom)) ;
359  }
360  else {
361  // Need to compute the derivative :
362  // Tensor for der
363  Tensor auxi_der (so.val_t->get_space(), val_res, type_ind, basis) ;
364  // Set the names of the indices :
365  auxi_der.set_name_affected() ;
366  auxi_der.set_name_ind(0, ind_der) ;
367  for (int i=1 ; i<val_res ; i++)
368  auxi_der.set_name_ind(i, so.der_t->get_name_ind()[i-1]) ;
369 
370  //Loop on the components :
371  Index pos_auxi_der(auxi_der) ;
372  do {
373  for (int i=0 ; i<val_res-1 ; i++)
374  pos_so.set(i) = pos_auxi_der(i+1) ;
375  auxi_der.set(pos_auxi_der).set_domain(num_dom) = (*so.der_t)(pos_so)(num_dom).der_abs(pos_auxi_der(0)+1) ;
376  }
377  while (pos_auxi_der.inc()) ;
378 
379  // Need for derivative :
380  Term_eq auxi (num_dom, auxi_val, auxi_der) ;
381 
382  // If derive contravariant : manipulate first indice :
383  if (type_der==CON)
384  manipulator->manipulate_ind (auxi, 0) ;
385 
386  if (!need_sum)
387  return auxi ;
388  else
389  return Term_eq (num_dom, auxi.val_t->do_summation_one_dom(num_dom),
390  auxi.der_t->do_summation_one_dom(num_dom)) ;
391  }
392 
393 }
394 
395 Term_eq Domain_shell::derive_flat_mtz (int type_der, char ind_der, const Term_eq& so, const Metric* manipulator) const {
396 
397  assert ((type_der==COV) || (type_der==CON)) ;
398  int val_res = so.val_t->get_valence() + 1 ;
399 
400  bool donames = (ind_der==' ') ? false : true ;
401  bool need_sum = false ;
402 
403  Array<int> type_ind (val_res) ;
404  type_ind.set(0) = COV ;
405  for (int i=1 ; i<val_res ; i++)
406  type_ind.set(i) = so.val_t->get_index_type(i-1) ;
407 
408  if (donames) {
409  if (so.val_t->get_valence()>0)
410  assert (so.val_t->is_name_affected()) ;
411  // Need for summation ?
412  for (int i=1 ; i<val_res ; i++)
413  if (ind_der== so.val_t->get_name_ind()[i-1])
414  need_sum = true ;
415  }
416 
417  // Tensor for val
418  Base_tensor basis (so.val_t->get_space(), MTZ_BASIS) ;
419  Tensor auxi_val (so.val_t->get_space(), val_res, type_ind, basis) ;
420 
421  if (donames) {
422  // Set the names of the indices :
423  auxi_val.set_name_affected() ;
424  auxi_val.set_name_ind(0, ind_der) ;
425  for (int i=1 ; i<val_res ; i++)
426  auxi_val.set_name_ind(i, so.val_t->get_name_ind()[i-1]) ;
427  }
428 
429  // Part derivative
430  //Loop on the components :
431  Index pos_auxi_bis(auxi_val) ;
432  Index pos_so_bis (*so.val_t) ;
433  do {
434  for (int i=0 ; i<val_res-1 ; i++)
435  pos_so_bis.set(i) = pos_auxi_bis(i+1) ;
436  switch (pos_auxi_bis(0)) {
437  case 0 :
438  // d/dr :
439  auxi_val.set(pos_auxi_bis).set_domain(num_dom) = (*so.val_t)(pos_so_bis)(num_dom).der_r() ;
440  break ;
441  case 1 :
442  // cost/r dtheta
443  auxi_val.set(pos_auxi_bis).set_domain(num_dom) = (*so.val_t)(pos_so_bis)(num_dom).der_var(2).div_r().mult_cos_theta() ;
444  break ;
445  case 2 :
446  // cost/r sint d/dphi
447  auxi_val.set(pos_auxi_bis).set_domain(num_dom) = (*so.val_t)(pos_so_bis)(num_dom).der_var(3).div_r().div_sin_theta().mult_cos_theta() ;
448  break ;
449  default :
450  cerr << "Bad indice in Domain_shell::derive_flat_mtz" << endl ;
451  abort() ;
452  }
453  }
454  while (pos_auxi_bis.inc()) ;
455 
456  // Loop indice summation on connection symbols
457  for (int ind_sum=0 ; ind_sum<val_res-1 ; ind_sum++) {
458 
459  //Loop on the components :
460  Index pos_auxi(auxi_val) ;
461  Index pos_so (*so.val_t) ;
462 
463  do {
464  for (int i=0 ; i<val_res-1 ; i++)
465  pos_so.set(i) = pos_auxi(i+1) ;
466  // Different cases of the derivative index :
467  switch (pos_auxi(0)) {
468  case 0 :
469  // Dr nothing
470  break ;
471  case 1 :
472  // Dtheta
473  // Different cases of the source index
474  switch (pos_auxi(ind_sum+1)) {
475  case 0 :
476  //Dtheta S_r
477  pos_so.set(ind_sum) = 1 ;
478  auxi_val.set(pos_auxi).set_domain(num_dom) -= (*so.val_t)(pos_so)(num_dom).div_r() ;
479  break ;
480  case 1 :
481  // Dtheta S_theta
482  pos_so.set(ind_sum) = 0 ;
483  auxi_val.set(pos_auxi).set_domain(num_dom) += (*so.val_t)(pos_so)(num_dom).div_r() ;
484  break ;
485  case 2 :
486  //Dtheta S_phi
487  break ;
488  default :
489  cerr << "Bad indice in Domain_shell::derive_flat_mtz" << endl ;
490  abort() ;
491  }
492  break ;
493  case 2 :
494  // Dphi
495  // Different cases of the source index
496  switch (pos_auxi(ind_sum+1)) {
497  case 0 :
498  //Dphi S_r
499  pos_so.set(ind_sum) = 2 ;
500  auxi_val.set(pos_auxi).set_domain(num_dom) -= (*so.val_t)(pos_so)(num_dom).div_r() ;
501  break ;
502  case 1 :
503  // Dphi S_theta
504  pos_so.set(ind_sum) = 2 ;
505  auxi_val.set(pos_auxi).set_domain(num_dom) -= (*so.val_t)(pos_so)(num_dom).div_r().div_sin_theta() ;
506  break ;
507  case 2 :
508  //Dphi S_phi
509  pos_so.set(ind_sum) = 0 ;
510  auxi_val.set(pos_auxi).set_domain(num_dom) += (*so.val_t)(pos_so)(num_dom).div_r() ;
511  pos_so.set(ind_sum) = 1 ;
512  auxi_val.set(pos_auxi).set_domain(num_dom) += (*so.val_t)(pos_so)(num_dom).div_r().div_sin_theta() ;
513  break ;
514  default :
515  cerr << "Bad indice in Domain_shell::derive_flat_mtz" << endl ;
516  abort() ;
517  }
518  break ;
519  default :
520  cerr << "Bad indice in Domain_shell::derive_flat_mtz" << endl ;
521  abort() ;
522  }
523  }
524  while (pos_auxi.inc()) ;
525  }
526 
527  if (so.der_t==0x0) {
528  // No need for derivative :
529  Term_eq auxi (num_dom, auxi_val) ;
530  // If derive contravariant : manipulate first indice :
531  if (type_der==CON)
532  manipulator->manipulate_ind (auxi, 0) ;
533 
534  if (!need_sum)
535  return auxi ;
536  else
537  return Term_eq (num_dom, auxi.val_t->do_summation_one_dom(num_dom)) ;
538  }
539  else {
540  // Need to compute the derivative :
541  // Tensor for der
542  Tensor auxi_der (so.val_t->get_space(), val_res, type_ind, basis) ;
543 
544  if (donames) {
545  // Set the names of the indices :
546  auxi_der.set_name_affected() ;
547  auxi_der.set_name_ind(0, ind_der) ;
548  for (int i=1 ; i<val_res ; i++)
549  auxi_der.set_name_ind(i, so.der_t->get_name_ind()[i-1]) ;
550  }
551 
552  // Part derivative
553  //Loop on the components :
554  Index pos_auxi_der_bis(auxi_der) ;
555  do {
556  for (int i=0 ; i<val_res-1 ; i++)
557  pos_so_bis.set(i) = pos_auxi_der_bis(i+1) ;
558  switch (pos_auxi_der_bis(0)) {
559  case 0 :
560  // d/dr :
561  auxi_der.set(pos_auxi_der_bis).set_domain(num_dom) = (*so.der_t)(pos_so_bis)(num_dom).der_r() ;
562  break ;
563  case 1 :
564  // cost/r dtheta
565  auxi_der.set(pos_auxi_der_bis).set_domain(num_dom) = (*so.der_t)(pos_so_bis)(num_dom).der_var(2).div_r().mult_cos_theta() ;
566  break ;
567  case 2 :
568  // cost/r sint d/dphi
569  auxi_der.set(pos_auxi_der_bis).set_domain(num_dom) = (*so.der_t)(pos_so_bis)(num_dom).der_var(3).div_r().div_sin_theta().mult_cos_theta() ;
570  break ;
571  default :
572  cerr << "Bad indice in Domain_shell::derive_flat_mtz" << endl ;
573  abort() ;
574  }
575  }
576  while (pos_auxi_der_bis.inc()) ;
577 
578  // Loop indice summation on connection symbols
579  for (int ind_sum=0 ; ind_sum<val_res-1 ; ind_sum++) {
580 
581  //Loop on the components :
582  Index pos_auxi_der(auxi_val) ;
583  Index pos_so (*so.val_t) ;
584 
585  do {
586  for (int i=0 ; i<val_res-1 ; i++)
587  pos_so.set(i) = pos_auxi_der(i+1) ;
588  // Different cases of the derivative index :
589  switch (pos_auxi_der(0)) {
590  case 0 :
591  // Dr nothing
592  break ;
593  case 1 :
594  // Dtheta
595  // Different cases of the source index
596  switch (pos_auxi_der(ind_sum+1)) {
597  case 0 :
598  //Dtheta S_r
599  pos_so.set(ind_sum) = 1 ;
600  auxi_der.set(pos_auxi_der).set_domain(num_dom) -= (*so.der_t)(pos_so)(num_dom).div_r() ;
601  break ;
602  case 1 :
603  // Dtheta S_theta
604  pos_so.set(ind_sum) = 0 ;
605  auxi_der.set(pos_auxi_der).set_domain(num_dom) += (*so.der_t)(pos_so)(num_dom).div_r() ;
606  break ;
607  case 2 :
608  //Dtheta S_phi
609  break ;
610  default :
611  cerr << "Bad indice in Domain_shell::derive_flat_mtz" << endl ;
612  abort() ;
613  }
614  break ;
615  case 2 :
616  // Dphi
617  // Different cases of the source index
618  switch (pos_auxi_der(ind_sum+1)) {
619  case 0 :
620  //Dphi S_r
621  pos_so.set(ind_sum) = 2 ;
622  auxi_der.set(pos_auxi_der).set_domain(num_dom) -= (*so.der_t)(pos_so)(num_dom).div_r() ;
623  break ;
624  case 1 :
625  // Dphi S_theta
626  pos_so.set(ind_sum) = 2 ;
627  auxi_der.set(pos_auxi_der).set_domain(num_dom) -= (*so.der_t)(pos_so)(num_dom).div_r().div_sin_theta() ;
628  break ;
629  case 2 :
630  //Dphi S_phi
631  pos_so.set(ind_sum) = 0 ;
632  auxi_der.set(pos_auxi_der).set_domain(num_dom) += (*so.der_t)(pos_so)(num_dom).div_r() ;
633  pos_so.set(ind_sum) = 1 ;
634  auxi_der.set(pos_auxi_der).set_domain(num_dom) += (*so.der_t)(pos_so)(num_dom).div_r().div_sin_theta() ;
635  break ;
636  default :
637  cerr << "Bad indice in Domain_shell::derive_flat_mtz" << endl ;
638  abort() ;
639  }
640  break ;
641  default :
642  cerr << "Bad indice in Domain_shell::derive_flat_mtz" << endl ;
643  abort() ;
644  }
645  }
646  while (pos_auxi_der.inc()) ;
647  }
648 
649 
650  // Need for derivative :
651  Term_eq auxi (num_dom, auxi_val, auxi_der) ;
652  // If derive contravariant : manipulate first indice :
653  if (type_der==CON)
654  manipulator->manipulate_ind (auxi, 0) ;
655  if (!need_sum)
656  return auxi ;
657  else
658  return Term_eq (num_dom, auxi.val_t->do_summation_one_dom(num_dom),
659  auxi.der_t->do_summation_one_dom(num_dom)) ;
660  }
661 }
662 
663 
664 }
Kadath::Array::set
reference set(const Index &pos)
Read/write of an element.
Definition: array.hpp:186
Kadath::Base_tensor
Describes the tensorial basis used by the various tensors.
Definition: base_tensor.hpp:49
Kadath::Domain_shell::div_r
virtual Val_domain div_r(const Val_domain &) const
Division by .
Definition: domain_shell_ope.cpp:212
Kadath::Domain_shell::derive_flat_cart
virtual Term_eq derive_flat_cart(int, char, const Term_eq &, const Metric *) const
Computes the flat derivative of a Term_eq, in Cartesian coordinates.
Definition: domain_shell_for_metric.cpp:299
Kadath::Domain_shell::der_r
virtual Val_domain der_r(const Val_domain &) const
Compute the radial derivative of a scalar field.
Definition: domain_shell_ope.cpp:218
Kadath::Domain_shell::derive_flat_spher
virtual Term_eq derive_flat_spher(int, char, const Term_eq &, const Metric *) const
Computes the flat derivative of a Term_eq, in spherical orthonormal coordinates.
Definition: domain_shell_for_metric.cpp:30
Kadath::Domain_shell::derive_flat_mtz
virtual Term_eq derive_flat_mtz(int, char, const Term_eq &, const Metric *) const
Computes the flat derivative of a Term_eq, in spherical coordinates where the constant radii sections...
Definition: domain_shell_for_metric.cpp:395
Kadath::Domain::num_dom
int num_dom
Number of the current domain (used by the Space)
Definition: space.hpp:63
Kadath::Index
Class that gives the position inside a multi-dimensional Array.
Definition: index.hpp:38
Kadath::Index::set
int & set(int i)
Read/write of the position in a given dimension.
Definition: index.hpp:72
Kadath::Index::inc
bool inc(int increm, int var=0)
Increments the position of the Index.
Definition: index.hpp:99
Kadath::Metric
Purely abstract class for metric handling.
Definition: metric.hpp:39
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::manipulate_ind
virtual void manipulate_ind(Term_eq &so, int ind) const
Uses the Metric to manipulate one of the index of a Term_eq (i.e.
Definition: metric.cpp:645
Kadath::Scalar::set_domain
Val_domain & set_domain(int)
Read/write of a particular Val_domain.
Definition: scalar.hpp:555
Kadath::Tensor
Tensor handling.
Definition: tensor.hpp:149
Kadath::Tensor::set_name_ind
void set_name_ind(int dd, char name)
Sets the name of one index ; the names must have been affected first.
Definition: tensor_impl.hpp:186
Kadath::Tensor::set_name_affected
void set_name_affected()
Affects the name of the indices.
Definition: tensor.hpp:435
Kadath::Tensor::set
Scalar & set(const Array< int > &ind)
Returns the value of a component (read/write version).
Definition: tensor_impl.hpp:91
Kadath::Tensor::get_name_ind
char const * get_name_ind() const
Definition: tensor.hpp:424
Kadath::Tensor::get_index_type
int get_index_type(int i) const
Gives the type (covariant or contravariant) of a given index.
Definition: tensor.hpp:526
Kadath::Tensor::get_valence
int get_valence() const
Returns the valence.
Definition: tensor.hpp:509
Kadath::Tensor::is_name_affected
bool is_name_affected() const
Check whether the names of the indices have been affected.
Definition: tensor.hpp:429
Kadath::Tensor::do_summation_one_dom
Tensor do_summation_one_dom(int dd) const
Does the inner contraction of the Tensor in a given domain.
Definition: tensor_math_one_dom.cpp:590
Kadath::Tensor::get_space
const Space & get_space() const
Returns the Space.
Definition: tensor.hpp:499
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::div_r
Val_domain div_r() const
Division by the radius.
Definition: val_domain.cpp:743
Kadath::Val_domain::div_sin_theta
Val_domain div_sin_theta() const
Division by .
Definition: val_domain_ope.cpp:50
Kadath::Val_domain::mult_cos_theta
Val_domain mult_cos_theta() const
Multiplication by .
Definition: val_domain_ope.cpp:36