LORENE
tensor_calculus_ext.C
1 /*
2  * Function external to class Tensor for tensor calculus
3  *
4  *
5  */
6 
7 /*
8  * Copyright (c) 2003-2004 Eric Gourgoulhon & Jerome Novak
9  *
10  * Copyright (c) 1999-2001 Philippe Grandclement (for preceding class Tenseur)
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License as published by
16  * the Free Software Foundation; either version 2 of the License, or
17  * (at your option) any later version.
18  *
19  * LORENE is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with LORENE; if not, write to the Free Software
26  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27  *
28  */
29 
30 
31 char tensor_calculus_ext_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/tensor_calculus_ext.C,v 1.14 2014/10/13 08:53:44 j_novak Exp $" ;
32 
33 /*
34  * $Id: tensor_calculus_ext.C,v 1.14 2014/10/13 08:53:44 j_novak Exp $
35  * $Log: tensor_calculus_ext.C,v $
36  * Revision 1.14 2014/10/13 08:53:44 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.13 2014/10/06 15:13:20 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.12 2008/12/05 08:44:02 j_novak
43  * New flag to control the "verbosity" of maxabs.
44  *
45  * Revision 1.11 2004/05/13 21:32:29 e_gourgoulhon
46  * Added functions central_value, max_all_domains,
47  * min_all_domains and maxabs_all_domains.
48  *
49  * Revision 1.10 2004/02/27 21:15:34 e_gourgoulhon
50  * Suppressed function contract_desal (since function contract has now
51  * the boolean argument "desaliasing").
52  *
53  * Revision 1.9 2004/02/19 22:11:00 e_gourgoulhon
54  * Added argument "comment" in functions max, min, etc...
55  *
56  * Revision 1.8 2004/02/18 18:51:04 e_gourgoulhon
57  * Tensorial product moved from file tensor_calculus.C, since
58  * it is not a method of class Tensor.
59  *
60  * Revision 1.7 2004/02/18 15:56:23 e_gourgoulhon
61  * -- Added function contract for double contraction.
62  * -- Efficiency improved in functions contract: better handling of variable
63  * "work"(work is now a reference on the relevant component of the result).
64  *
65  * Revision 1.6 2004/01/23 08:00:16 e_gourgoulhon
66  * Minor modifs. in output of methods min, max, maxabs, diffrel to
67  * better handle the display in the scalar case.
68  *
69  * Revision 1.5 2004/01/15 10:59:53 f_limousin
70  * Added method contract_desal for the contraction of two tensors with desaliasing
71  *
72  * Revision 1.4 2004/01/14 11:38:32 f_limousin
73  * Added method contract for one tensor
74  *
75  * Revision 1.3 2003/11/05 15:29:36 e_gourgoulhon
76  * Added declaration of externa functions max, min, maxabs,
77  * diffrel and diffrelmax.
78  *
79  * Revision 1.2 2003/10/11 16:47:10 e_gourgoulhon
80  * Suppressed the call to Ibtl::set_etat_qcq() after the construction
81  * of the Itbl's, thanks to the new property of the Itbl class.
82  *
83  * Revision 1.1 2003/10/06 15:13:38 e_gourgoulhon
84  * Tensor contraction.
85  *
86  *
87  * $Header: /cvsroot/Lorene/C++/Source/Tensor/tensor_calculus_ext.C,v 1.14 2014/10/13 08:53:44 j_novak Exp $
88  *
89  */
90 
91 // Headers C
92 #include <cstdlib>
93 #include <cassert>
94 #include <cmath>
95 
96 // Headers Lorene
97 #include "tensor.h"
98 
99 
100 
101  //-----------------------//
102  // Tensorial product //
103  //-----------------------//
104 
105 
106 namespace Lorene {
107 Tensor operator*(const Tensor& t1, const Tensor& t2) {
108 
109  assert (t1.mp == t2.mp) ;
110 
111  int val_res = t1.valence + t2.valence ;
112 
113  Itbl tipe (val_res) ;
114 
115  for (int i=0 ; i<t1.valence ; i++)
116  tipe.set(i) = t1.type_indice(i) ;
117  for (int i=0 ; i<t2.valence ; i++)
118  tipe.set(i+t1.valence) = t2.type_indice(i) ;
119 
120 
121  if ( (t1.valence != 0) && (t2.valence != 0) ) {
122  assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
123  }
124 
125  const Base_vect* triad_res ;
126  if (t1.valence != 0) {
127  triad_res = t1.get_triad() ;
128  }
129  else{
130  triad_res = t2.get_triad() ;
131  }
132 
133  Tensor res(*t1.mp, val_res, tipe, triad_res) ;
134 
135  Itbl jeux_indice_t1 (t1.valence) ;
136  Itbl jeux_indice_t2 (t2.valence) ;
137 
138  for (int i=0 ; i<res.n_comp ; i++) {
139  Itbl jeux_indice_res(res.indices(i)) ;
140  for (int j=0 ; j<t1.valence ; j++)
141  jeux_indice_t1.set(j) = jeux_indice_res(j) ;
142  for (int j=0 ; j<t2.valence ; j++)
143  jeux_indice_t2.set(j) = jeux_indice_res(j+t1.valence) ;
144 
145  res.set(jeux_indice_res) = t1(jeux_indice_t1)*t2(jeux_indice_t2) ;
146  }
147 
148  return res ;
149 }
150 
151 
152 
153 
154  //------------------//
155  // Contraction //
156  //------------------//
157 
158 // Contraction on one index
159 // ------------------------
160 Tensor contract(const Tensor& t1, int ind1, const Tensor& t2, int ind2,
161  bool desaliasing) {
162 
163  int val1 = t1.get_valence() ;
164  int val2 = t2.get_valence() ;
165 
166  // Verifs :
167  assert((ind1>=0) && (ind1<val1)) ;
168  assert((ind2>=0) && (ind2<val2)) ;
169  assert(t1.get_mp() == t2.get_mp()) ;
170 
171  // Contraction possible ?
172  if ( (val1 != 0) && (val2 != 0) ) {
173  assert ( *(t1.get_triad()) == *(t2.get_triad()) ) ;
174  }
175  assert (t1.get_index_type(ind1) != t2.get_index_type(ind2)) ;
176 
177  int val_res = val1 + val2 - 2;
178 
179  Itbl tipe(val_res) ;
180 
181  for (int i=0 ; i<ind1 ; i++)
182  tipe.set(i) = t1.get_index_type(i) ;
183  for (int i=ind1 ; i<val1-1 ; i++)
184  tipe.set(i) = t1.get_index_type(i+1) ;
185  for (int i=val1-1 ; i<val1+ind2-1 ; i++)
186  tipe.set(i) = t2.get_index_type(i-val1+1) ;
187  for (int i = val1+ind2-1 ; i<val_res ; i++)
188  tipe.set(i) = t2.get_index_type(i-val1+2) ;
189 
190  const Base_vect* triad_res = (val_res == 0) ? 0x0 : t1.get_triad() ;
191 
192  Tensor res(t1.get_mp(), val_res, tipe, triad_res) ;
193 
194  // Boucle sur les composantes de res :
195 
196  Itbl jeux_indice_t1(val1) ;
197  Itbl jeux_indice_t2(val2) ;
198 
199  for (int i=0 ; i<res.get_n_comp() ; i++) {
200 
201  Itbl jeux_indice_res(res.indices(i)) ;
202 
203  for (int j=0 ; j<ind1 ; j++)
204  jeux_indice_t1.set(j) = jeux_indice_res(j) ;
205 
206  for (int j=ind1+1 ; j<val1 ; j++)
207  jeux_indice_t1.set(j) = jeux_indice_res(j-1) ;
208 
209  for (int j=0 ; j<ind2 ; j++)
210  jeux_indice_t2.set(j) = jeux_indice_res(val1+j-1) ;
211 
212  for (int j=ind2+1 ; j<val2 ; j++)
213  jeux_indice_t2.set(j) = jeux_indice_res(val1+j-2) ;
214 
215  Scalar& work = res.set(jeux_indice_res) ;
216  work.set_etat_zero() ;
217 
218  for (int j=1 ; j<=3 ; j++) {
219  jeux_indice_t1.set(ind1) = j ;
220  jeux_indice_t2.set(ind2) = j ;
221  if (desaliasing) {
222  work += t1(jeux_indice_t1) % t2(jeux_indice_t2) ;
223  }
224  else {
225  work += t1(jeux_indice_t1) * t2(jeux_indice_t2) ;
226  }
227  }
228 
229  }
230 
231  return res ;
232 }
233 
234 
235 
236 // Contraction on two indices
237 // --------------------------
238 Tensor contract(const Tensor& t1, int i1, int j1,
239  const Tensor& t2, int i2, int j2,
240  bool desaliasing) {
241 
242  int val1 = t1.get_valence() ;
243  int val2 = t2.get_valence() ;
244 
245  // Verifs :
246  assert( val1 >= 2 ) ;
247  assert( val2 >= 2 ) ;
248  assert( (0<=i1) && (i1<j1) && (j1<val1) ) ;
249  assert( (0<=i2) && (i2<j2) && (j2<val2) ) ;
250  assert( t1.get_mp() == t2.get_mp() ) ;
251 
252  // Contraction possible ?
253  assert( *(t1.get_triad()) == *(t2.get_triad()) ) ;
254  assert( t1.get_index_type(i1) != t2.get_index_type(i2) ) ;
255  assert( t1.get_index_type(j1) != t2.get_index_type(j2) ) ;
256 
257  int val_res = val1 + val2 - 4 ;
258 
259  Itbl tipe(val_res) ;
260 
261  for (int i=0 ; i<i1 ; i++)
262  tipe.set(i) = t1.get_index_type(i) ;
263 
264  for (int i=i1 ; i<j1-1 ; i++)
265  tipe.set(i) = t1.get_index_type(i+1) ;
266 
267  for (int i=j1-1 ; i<val1-2 ; i++)
268  tipe.set(i) = t1.get_index_type(i+2) ;
269 
270  for (int i=val1-2 ; i<val1-2+i2 ; i++)
271  tipe.set(i) = t2.get_index_type(i-val1+2) ;
272 
273  for (int i=val1-2+i2 ; i<val1+j2-3 ; i++)
274  tipe.set(i) = t2.get_index_type(i-val1+3) ;
275 
276  for (int i=val1+j2-3 ; i<val_res ; i++)
277  tipe.set(i) = t2.get_index_type(i-val1+4) ;
278 
279  const Base_vect* triad_res = (val_res == 0) ? 0x0 : t1.get_triad() ;
280 
281  Tensor res(t1.get_mp(), val_res, tipe, triad_res) ;
282 
283  // Boucle sur les composantes de res :
284 
285  Itbl jeux_indice_t1(val1) ;
286  Itbl jeux_indice_t2(val2) ;
287 
288  for (int ic=0 ; ic<res.get_n_comp() ; ic++) {
289 
290  Itbl jeux_indice_res(res.indices(ic)) ;
291 
292  for (int k=0 ; k<i1 ; k++)
293  jeux_indice_t1.set(k) = jeux_indice_res(k) ;
294 
295  for (int k=i1+1 ; k<j1 ; k++)
296  jeux_indice_t1.set(k) = jeux_indice_res(k-1) ;
297 
298  for (int k=j1+1 ; k<val1 ; k++)
299  jeux_indice_t1.set(k) = jeux_indice_res(k-2) ;
300 
301  for (int k=0 ; k<i2 ; k++)
302  jeux_indice_t2.set(k) = jeux_indice_res(val1+k-2) ;
303 
304  for (int k=i2+1 ; k<j2 ; k++)
305  jeux_indice_t2.set(k) = jeux_indice_res(val1+k-3) ;
306 
307  for (int k=j2+1 ; k<val2 ; k++)
308  jeux_indice_t2.set(k) = jeux_indice_res(val1+k-4) ;
309 
310  Scalar& work = res.set(jeux_indice_res) ;
311  work.set_etat_zero() ;
312 
313  for (int i=1 ; i<=3 ; i++) {
314 
315  jeux_indice_t1.set(i1) = i ;
316  jeux_indice_t2.set(i2) = i ;
317 
318  for (int j=1 ; j<=3 ; j++) {
319 
320  jeux_indice_t1.set(j1) = j ;
321  jeux_indice_t2.set(j2) = j ;
322 
323  if (desaliasing) {
324  work += t1(jeux_indice_t1) % t2(jeux_indice_t2) ;
325  }
326  else {
327  work += t1(jeux_indice_t1) * t2(jeux_indice_t2) ;
328  }
329  }
330  }
331 
332  }
333 
334  return res ;
335 }
336 
337 
338 
339 
340 Tensor contract(const Tensor& source, int ind_1, int ind_2) {
341 
342  int val = source.get_valence() ;
343 
344  // Les verifications :
345  assert ((ind_1 >= 0) && (ind_1 < val)) ;
346  assert ((ind_2 >= 0) && (ind_2 < val)) ;
347  assert (ind_1 != ind_2) ;
348  assert (source.get_index_type(ind_1) != source.get_index_type(ind_2)) ;
349 
350  // On veut ind_1 < ind_2 :
351  if (ind_1 > ind_2) {
352  int auxi = ind_2 ;
353  ind_2 = ind_1 ;
354  ind_1 = auxi ;
355  }
356 
357  // On construit le resultat :
358  int val_res = val - 2 ;
359 
360  Itbl tipe(val_res) ;
361 
362  for (int i=0 ; i<ind_1 ; i++)
363  tipe.set(i) = source.get_index_type(i) ;
364  for (int i=ind_1 ; i<ind_2-1 ; i++)
365  tipe.set(i) = source.get_index_type(i+1) ;
366  for (int i = ind_2-1 ; i<val_res ; i++)
367  tipe.set(i) = source.get_index_type(i+2) ;
368 
369  Tensor res(source.get_mp(), val_res, tipe, source.get_triad()) ;
370 
371  // Boucle sur les composantes de res :
372 
373  Itbl jeux_indice_source(val) ;
374 
375  for (int i=0 ; i<res.get_n_comp() ; i++) {
376 
377  Itbl jeux_indice_res(res.indices(i)) ;
378 
379  for (int j=0 ; j<ind_1 ; j++)
380  jeux_indice_source.set(j) = jeux_indice_res(j) ;
381  for (int j=ind_1+1 ; j<ind_2 ; j++)
382  jeux_indice_source.set(j) = jeux_indice_res(j-1) ;
383  for (int j=ind_2+1 ; j<val ; j++)
384  jeux_indice_source.set(j) = jeux_indice_res(j-2) ;
385 
386  Scalar& work = res.set(jeux_indice_res) ;
387  work.set_etat_zero() ;
388 
389  for (int j=1 ; j<=3 ; j++) {
390  jeux_indice_source.set(ind_1) = j ;
391  jeux_indice_source.set(ind_2) = j ;
392  work += source(jeux_indice_source) ;
393  }
394 
395  }
396 
397  return res ;
398 }
399 
400 
401 
402 
403  //------------------//
404  // diffrel //
405  //------------------//
406 
407 
408 Tbl diffrel(const Tensor& aa, const Tensor& bb, const char* comment,
409  ostream& ost) {
410 
411  if (comment != 0x0) ost << comment << " : " << endl ;
412 
413  int val = aa.get_valence() ;
414 
415  assert(bb.get_valence() == val) ;
416 
417  int n_comp_a = aa.get_n_comp() ;
418  int n_comp_b = bb.get_n_comp() ;
419 
420  const Tensor* tmax ;
421  int n_comp_max ;
422  if (n_comp_a >= n_comp_b) {
423  n_comp_max = n_comp_a ;
424  tmax = &aa ;
425  }
426  else {
427  n_comp_max = n_comp_b ;
428  tmax = &bb ;
429  }
430 
431  int nz = aa.get_mp().get_mg()->get_nzone() ;
432  Tbl resu(n_comp_max, nz) ;
433  resu.set_etat_qcq() ;
434 
435  Itbl idx(val) ;
436 
437  for (int ic=0; ic<n_comp_max; ic++) {
438  idx = tmax->indices(ic) ;
439  Tbl diff = diffrel( aa(idx), bb(idx) ) ;
440 
441  if (n_comp_max > 1) ost << " Comp." ;
442  for (int j=0 ; j<val ; j++) {
443  ost << " " << idx(j) ;
444  }
445  if (n_comp_max > 1) ost << " : " ;
446  else ost << " " ;
447  for (int l=0; l<nz; l++) {
448  ost << " " << diff(l) ;
449  resu.set(ic, l) = diff(l) ;
450  }
451  ost << "\n" ;
452 
453  }
454 
455  return resu ;
456 }
457 
458 
459  //--------------------//
460  // diffrelmax //
461  //--------------------//
462 
463 
464 Tbl diffrelmax(const Tensor& aa, const Tensor& bb, const char* comment,
465  ostream& ost) {
466 
467  if (comment != 0x0) ost << comment << " : " << endl ;
468 
469  int val = aa.get_valence() ;
470 
471  assert(bb.get_valence() == val) ;
472 
473  int n_comp_a = aa.get_n_comp() ;
474  int n_comp_b = bb.get_n_comp() ;
475 
476  const Tensor* tmax ;
477  int n_comp_max ;
478  if (n_comp_a >= n_comp_b) {
479  n_comp_max = n_comp_a ;
480  tmax = &aa ;
481  }
482  else {
483  n_comp_max = n_comp_b ;
484  tmax = &bb ;
485  }
486 
487  int nz = aa.get_mp().get_mg()->get_nzone() ;
488  Tbl resu(n_comp_max, nz) ;
489  resu.set_etat_qcq() ;
490 
491  Itbl idx(val) ;
492 
493  for (int ic=0; ic<n_comp_max; ic++) {
494  idx = tmax->indices(ic) ;
495  Tbl diff = diffrelmax( aa(idx), bb(idx) ) ;
496 
497  if (n_comp_max > 1) ost << " Comp." ;
498  for (int j=0 ; j<val ; j++) {
499  ost << " " << idx(j) ;
500  }
501  if (n_comp_max > 1) ost << " : " ;
502  else ost << " " ;
503  for (int l=0; l<nz; l++) {
504  ost << " " << diff(l) ;
505  resu.set(ic, l) = diff(l) ;
506  }
507  ost << "\n" ;
508 
509  }
510 
511  return resu ;
512 }
513 
514 
515 
516  //----------------//
517  // max //
518  //----------------//
519 
520 
521 Tbl max(const Tensor& aa, const char* comment, ostream& ost) {
522 
523  if (comment != 0x0) ost << comment << " : " << endl ;
524 
525  int val = aa.get_valence() ;
526 
527  int n_comp = aa.get_n_comp() ;
528 
529  int nz = aa.get_mp().get_mg()->get_nzone() ;
530  Tbl resu(n_comp, nz) ;
531  resu.set_etat_qcq() ;
532 
533  Itbl idx(val) ;
534 
535  for (int ic=0; ic<n_comp; ic++) {
536 
537  idx = aa.indices(ic) ;
538  Tbl diff = max( aa(idx) ) ;
539 
540  if (val > 0) ost << " Comp." ;
541  for (int j=0 ; j<val ; j++) {
542  ost << " " << idx(j) ;
543  }
544  if (val > 0) ost << " : " ;
545  else ost << " " ;
546  for (int l=0; l<nz; l++) {
547  ost << " " << diff(l) ;
548  resu.set(ic, l) = diff(l) ;
549  }
550  ost << "\n" ;
551 
552  }
553 
554  return resu ;
555 }
556 
557 
558 
559  //----------------//
560  // min //
561  //----------------//
562 
563 
564 Tbl min(const Tensor& aa, const char* comment, ostream& ost) {
565 
566  if (comment != 0x0) ost << comment << " : " << endl ;
567 
568  int val = aa.get_valence() ;
569 
570  int n_comp = aa.get_n_comp() ;
571 
572  int nz = aa.get_mp().get_mg()->get_nzone() ;
573  Tbl resu(n_comp, nz) ;
574  resu.set_etat_qcq() ;
575 
576  Itbl idx(val) ;
577 
578  for (int ic=0; ic<n_comp; ic++) {
579 
580  idx = aa.indices(ic) ;
581  Tbl diff = min( aa(idx) ) ;
582 
583  if (val > 0) ost << " Comp." ;
584  for (int j=0 ; j<val ; j++) {
585  ost << " " << idx(j) ;
586  }
587  if (val > 0) ost << " : " ;
588  else ost << " " ;
589  for (int l=0; l<nz; l++) {
590  ost << " " << diff(l) ;
591  resu.set(ic, l) = diff(l) ;
592  }
593  ost << "\n" ;
594 
595  }
596 
597  return resu ;
598 }
599 
600 
601  //--------------------//
602  // maxabs //
603  //--------------------//
604 
605 
606 Tbl maxabs(const Tensor& aa, const char* comment, ostream& ost, bool verb) {
607 
608  if (comment != 0x0) ost << comment << " : " << endl ;
609 
610  int val = aa.get_valence() ;
611 
612  int n_comp = aa.get_n_comp() ;
613 
614  int nz = aa.get_mp().get_mg()->get_nzone() ;
615  Tbl resu(n_comp, nz) ;
616  resu.set_etat_qcq() ;
617 
618  Itbl idx(val) ;
619 
620  for (int ic=0; ic<n_comp; ic++) {
621 
622  idx = aa.indices(ic) ;
623  Tbl diff = max( abs( aa(idx) ) ) ;
624 
625  if (verb) {
626  if (val > 0) ost << " Comp." ;
627  for (int j=0 ; j<val ; j++) {
628  ost << " " << idx(j) ;
629  }
630  if (val > 0 ) ost << " : " ;
631  else ost << " " ;
632  }
633  for (int l=0; l<nz; l++) {
634  if (verb) ost << " " << diff(l) ;
635  resu.set(ic, l) = diff(l) ;
636  }
637  if (verb) ost << "\n" ;
638 
639  }
640 
641  return resu ;
642 }
643 
644 
645  //-------------------//
646  // Central value //
647  //-------------------//
648 
649 Tbl central_value(const Tensor& aa, const char* comment, ostream& ost) {
650 
651  if (comment != 0x0) ost << comment << " : " << endl ;
652 
653  int val = aa.get_valence() ;
654  int n_comp = aa.get_n_comp() ;
655 
656  Tbl resu(n_comp) ;
657  resu.set_etat_qcq() ;
658 
659  Itbl idx(val) ;
660 
661  for (int ic=0; ic<n_comp; ic++) {
662 
663  idx = aa.indices(ic) ;
664  double aa_c = aa(idx).val_grid_point(0,0,0,0) ;
665  resu.set(ic) = aa_c ;
666 
667  if ( comment != 0x0 ) {
668  if ( val > 0 ) ost << " Comp." ;
669  for (int j=0 ; j<val ; j++) {
670  ost << " " << idx(j) ;
671  }
672  if (val > 0 ) ost << " : " ;
673  else ost << " " ;
674  ost << aa_c << endl ;
675  }
676 
677  }
678 
679  return resu ;
680 }
681 
682 
683  //-------------------//
684  // max_all_domains //
685  //-------------------//
686 
687 Tbl max_all_domains(const Tensor& aa, int l_excluded, const char* comment,
688  ostream& ost) {
689 
690  if (comment != 0x0) ost << comment << " : " << endl ;
691 
692  Tbl max_dom = max(aa) ;
693 
694  int val = aa.get_valence() ;
695  int n_comp = aa.get_n_comp() ;
696  int nz = aa.get_mp().get_mg()->get_nzone() ;
697 
698  Tbl resu(n_comp) ;
699  resu.set_etat_qcq() ;
700 
701  Itbl idx(val) ;
702 
703  for (int ic=0; ic<n_comp; ic++) {
704 
705  double x0 ;
706  if (l_excluded != 0) x0 = max_dom(ic, 0) ;
707  else x0 = max_dom(ic, 1) ;
708  for (int l=0; l<nz; l++) {
709  if (l == l_excluded) continue ;
710  double x = max_dom(ic,l) ;
711  if (x > x0) x0 = x ;
712  }
713 
714  resu.set(ic) = x0 ;
715 
716  if ( comment != 0x0 ) {
717  if ( val > 0 ) ost << " Comp." ;
718  idx = aa.indices(ic) ;
719  for (int j=0 ; j<val ; j++) {
720  ost << " " << idx(j) ;
721  }
722  if (val > 0 ) ost << " : " ;
723  else ost << " " ;
724  ost << x0 << endl ;
725  }
726 
727  }
728 
729  return resu ;
730 
731 }
732 
733  //-------------------//
734  // min_all_domains //
735  //-------------------//
736 
737 Tbl min_all_domains(const Tensor& aa, int l_excluded, const char* comment,
738  ostream& ost) {
739 
740  if (comment != 0x0) ost << comment << " : " << endl ;
741 
742  Tbl min_dom = min(aa) ;
743 
744  int val = aa.get_valence() ;
745  int n_comp = aa.get_n_comp() ;
746  int nz = aa.get_mp().get_mg()->get_nzone() ;
747 
748  Tbl resu(n_comp) ;
749  resu.set_etat_qcq() ;
750 
751  Itbl idx(val) ;
752 
753  for (int ic=0; ic<n_comp; ic++) {
754 
755  double x0 ;
756  if (l_excluded != 0) x0 = min_dom(ic, 0) ;
757  else x0 = min_dom(ic, 1) ;
758  for (int l=0; l<nz; l++) {
759  if (l == l_excluded) continue ;
760  double x = min_dom(ic,l) ;
761  if (x < x0) x0 = x ;
762  }
763 
764  resu.set(ic) = x0 ;
765 
766  if ( comment != 0x0 ) {
767  if ( val > 0 ) ost << " Comp." ;
768  idx = aa.indices(ic) ;
769  for (int j=0 ; j<val ; j++) {
770  ost << " " << idx(j) ;
771  }
772  if (val > 0 ) ost << " : " ;
773  else ost << " " ;
774  ost << x0 << endl ;
775  }
776 
777  }
778 
779  return resu ;
780 
781 }
782 
783 
784  //----------------------//
785  // maxabs_all_domains //
786  //----------------------//
787 
788 Tbl maxabs_all_domains(const Tensor& aa, int l_excluded, const char* comment,
789  ostream& ost, bool verb) {
790 
791  if (comment != 0x0) ost << comment << " : " << endl ;
792 
793  Tbl maxabs_dom = maxabs(aa, 0x0, ost, verb) ;
794 
795  int val = aa.get_valence() ;
796  int n_comp = aa.get_n_comp() ;
797  int nz = aa.get_mp().get_mg()->get_nzone() ;
798 
799  Tbl resu(n_comp) ;
800  resu.set_etat_qcq() ;
801 
802  Itbl idx(val) ;
803 
804  for (int ic=0; ic<n_comp; ic++) {
805 
806  double x0 ;
807  if (l_excluded != 0) x0 = maxabs_dom(ic, 0) ;
808  else x0 = maxabs_dom(ic, 1) ;
809  for (int l=0; l<nz; l++) {
810  if (l == l_excluded) continue ;
811  double x = maxabs_dom(ic,l) ;
812  if (x > x0) x0 = x ;
813  }
814 
815  resu.set(ic) = x0 ;
816 
817  if ( comment != 0x0 ) {
818  if ( val > 0 ) ost << " Comp." ;
819  idx = aa.indices(ic) ;
820  for (int j=0 ; j<val ; j++) {
821  ost << " " << idx(j) ;
822  }
823  if (val > 0 ) ost << " : " ;
824  else ost << " " ;
825  ost << x0 << endl ;
826  }
827 
828  }
829 
830  return resu ;
831 
832 }
833 
834 
835 
836 }
Tbl maxabs(const Tensor &aa, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maxima in each domain of the absolute values of the tensor components.
Tbl min_all_domains(const Tensor &aa, int l_excluded=-1, const char *comment=0x0, ostream &ost=cout)
Minimum value of each component of a tensor over all the domains.
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:312
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:324
Tbl central_value(const Tensor &aa, const char *comment=0x0, ostream &ost=cout)
Central value of each component of a tensor.
Lorene prototypes.
Definition: app_hor.h:64
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Base_val operator*(const Base_val &, const Base_val &)
This operator is used when calling multiplication or division of Valeur .
int get_n_comp() const
Returns the number of stored components.
Definition: tensor.h:872
Basic integer array class.
Definition: itbl.h:122
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition: cmp_math.C:458
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition: base_vect.h:105
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Itbl type_indice
1D array of integers (class Itbl ) of size valence containing the type of each index: COV for a covar...
Definition: tensor.h:310
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:866
Tbl maxabs_all_domains(const Tensor &aa, int l_excluded=-1, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maximum of the absolute value of each component of a tensor over all the domains. ...
Tbl max_all_domains(const Tensor &aa, int l_excluded=-1, const char *comment=0x0, ostream &ost=cout)
Maximum value of each component of a tensor over all the domains.
int get_index_type(int i) const
Gives the type (covariant or contravariant) of the index number i .
Definition: tensor.h:886
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
Tensor handling.
Definition: tensor.h:288
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
int get_valence() const
Returns the valence.
Definition: tensor.h:869
int valence
Valence of the tensor (0 = scalar, 1 = vector, etc...)
Definition: tensor.h:298
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:410
Basic array class.
Definition: tbl.h:161
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:295
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:861
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:539
virtual Itbl indices(int pos) const
Returns the indices of a component given by its position in the array cmp .
Definition: tensor.C:539