ginac/tensor.cpp

   1 /** @file tensor.cpp
   2  *
   3  *  Implementation of GiNaC's special tensors. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
   7  *
   8  *  This program is free software; you can redistribute it and/or modify
   9  *  it under the terms of the GNU General Public License as published by
  10  *  the Free Software Foundation; either version 2 of the License, or
  11  *  (at your option) any later version.
  12  *
  13  *  This program is distributed in the hope that it will be useful,
  14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  *  GNU General Public License for more details.
  17  *
  18  *  You should have received a copy of the GNU General Public License
  19  *  along with this program; if not, write to the Free Software
  20  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  21  */
  22
  23 #include <stdexcept>
  24 #include <vector>
  25
  26 #include "tensor.h"
  27 #include "idx.h"
  28 #include "indexed.h"
  29 #include "relational.h"
  30 #include "lst.h"
  31 #include "numeric.h"
  32 #include "print.h"
  33 #include "archive.h"
  34 #include "utils.h"
  35 #include "debugmsg.h"
  36
  37 namespace GiNaC {
  38
  39 GINAC_IMPLEMENT_REGISTERED_CLASS(tensor, basic)
  40 GINAC_IMPLEMENT_REGISTERED_CLASS(tensdelta, tensor)
  41 GINAC_IMPLEMENT_REGISTERED_CLASS(tensmetric, tensor)
  42 GINAC_IMPLEMENT_REGISTERED_CLASS(minkmetric, tensmetric)
  43 GINAC_IMPLEMENT_REGISTERED_CLASS(tensepsilon, tensor)
  44
  45 //////////
  46 // default constructor, destructor, copy constructor assignment operator and helpers
  47 //////////
  48
  49 tensor::tensor(unsigned ti) : inherited(ti)
  50 {
  51         debugmsg("tensor constructor from unsigned", LOGLEVEL_CONSTRUCT); \
  52 }
  53
  54 DEFAULT_CTORS(tensor)
  55 DEFAULT_CTORS(tensdelta)
  56 DEFAULT_CTORS(tensmetric)
  57 DEFAULT_DESTROY(minkmetric)
  58 DEFAULT_DESTROY(tensepsilon)
  59
  60 minkmetric::minkmetric() : pos_sig(false)
  61 {
  62         debugmsg("minkmetric default constructor", LOGLEVEL_CONSTRUCT);
  63         tinfo_key = TINFO_minkmetric;
  64 }
  65
  66 minkmetric::minkmetric(bool ps) : pos_sig(ps)
  67 {
  68         debugmsg("minkmetric constructor from bool", LOGLEVEL_CONSTRUCT);
  69         tinfo_key = TINFO_minkmetric;
  70 }
  71
  72 void minkmetric::copy(const minkmetric & other)
  73 {
  74         inherited::copy(other);
  75         pos_sig = other.pos_sig;
  76 }
  77
  78 tensepsilon::tensepsilon() : minkowski(false), pos_sig(false)
  79 {
  80         debugmsg("tensepsilon default constructor", LOGLEVEL_CONSTRUCT);
  81         tinfo_key = TINFO_tensepsilon;
  82 }
  83
  84 tensepsilon::tensepsilon(bool mink, bool ps) : minkowski(mink), pos_sig(ps)
  85 {
  86         debugmsg("tensepsilon constructor from bool,bool", LOGLEVEL_CONSTRUCT);
  87         tinfo_key = TINFO_tensepsilon;
  88 }
  89
  90 void tensepsilon::copy(const tensepsilon & other)
  91 {
  92         inherited::copy(other);
  93         minkowski = other.minkowski;
  94         pos_sig = other.pos_sig;
  95 }
  96
  97 //////////
  98 // archiving
  99 //////////
 100
 101 DEFAULT_ARCHIVING(tensor)
 102 DEFAULT_ARCHIVING(tensdelta)
 103 DEFAULT_ARCHIVING(tensmetric)
 104 DEFAULT_UNARCHIVE(minkmetric)
 105 DEFAULT_UNARCHIVE(tensepsilon)
 106
 107 minkmetric::minkmetric(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
 108 {
 109         debugmsg("minkmetric constructor from archive_node", LOGLEVEL_CONSTRUCT);
 110         n.find_bool("pos_sig", pos_sig);
 111 }
 112
 113 void minkmetric::archive(archive_node &n) const
 114 {
 115         inherited::archive(n);
 116         n.add_bool("pos_sig", pos_sig);
 117 }
 118
 119 tensepsilon::tensepsilon(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
 120 {
 121         debugmsg("tensepsilon constructor from archive_node", LOGLEVEL_CONSTRUCT);
 122         n.find_bool("minkowski", minkowski);
 123         n.find_bool("pos_sig", pos_sig);
 124 }
 125
 126 void tensepsilon::archive(archive_node &n) const
 127 {
 128         inherited::archive(n);
 129         n.add_bool("minkowski", minkowski);
 130         n.add_bool("pos_sig", pos_sig);
 131 }
 132
 133 //////////
 134 // functions overriding virtual functions from bases classes
 135 //////////
 136
 137 DEFAULT_COMPARE(tensor)
 138 DEFAULT_COMPARE(tensdelta)
 139 DEFAULT_COMPARE(tensmetric)
 140
 141 int minkmetric::compare_same_type(const basic & other) const
 142 {
 143         GINAC_ASSERT(is_of_type(other, minkmetric));
 144         const minkmetric &o = static_cast<const minkmetric &>(other);
 145
 146         if (pos_sig != o.pos_sig)
 147                 return pos_sig ? -1 : 1;
 148         else
 149                 return inherited::compare_same_type(other);
 150 }
 151
 152 int tensepsilon::compare_same_type(const basic & other) const
 153 {
 154         GINAC_ASSERT(is_of_type(other, tensepsilon));
 155         const tensepsilon &o = static_cast<const tensepsilon &>(other);
 156
 157         if (minkowski != o.minkowski)
 158                 return minkowski ? -1 : 1;
 159         else if (pos_sig != o.pos_sig)
 160                 return pos_sig ? -1 : 1;
 161         else
 162                 return inherited::compare_same_type(other);
 163 }
 164
 165 DEFAULT_PRINT(tensdelta, "delta")
 166 DEFAULT_PRINT(tensmetric, "g")
 167 DEFAULT_PRINT(minkmetric, "eta")
 168 DEFAULT_PRINT(tensepsilon, "eps")
 169
 170 /** Automatic symbolic evaluation of an indexed delta tensor. */
 171 ex tensdelta::eval_indexed(const basic & i) const
 172 {
 173         GINAC_ASSERT(is_of_type(i, indexed));
 174         GINAC_ASSERT(i.nops() == 3);
 175         GINAC_ASSERT(is_ex_of_type(i.op(0), tensdelta));
 176
 177         const idx & i1 = ex_to_idx(i.op(1));
 178         const idx & i2 = ex_to_idx(i.op(2));
 179
 180         // Trace of delta tensor is the dimension of the space
 181         if (is_dummy_pair(i1, i2))
 182                 return i1.get_dim();
 183
 184         // Numeric evaluation
 185         if (static_cast<const indexed &>(i).all_index_values_are(info_flags::integer)) {
 186                 int n1 = ex_to_numeric(i1.get_value()).to_int(), n2 = ex_to_numeric(i2.get_value()).to_int();
 187                 if (n1 == n2)
 188                         return _ex1();
 189                 else
 190                         return _ex0();
 191         }
 192
 193         // No further simplifications
 194         return i.hold();
 195 }
 196
 197 /** Automatic symbolic evaluation of an indexed metric tensor. */
 198 ex tensmetric::eval_indexed(const basic & i) const
 199 {
 200         GINAC_ASSERT(is_of_type(i, indexed));
 201         GINAC_ASSERT(i.nops() == 3);
 202         GINAC_ASSERT(is_ex_of_type(i.op(0), tensmetric));
 203         GINAC_ASSERT(is_ex_of_type(i.op(1), varidx));
 204         GINAC_ASSERT(is_ex_of_type(i.op(2), varidx));
 205
 206         const varidx & i1 = ex_to_varidx(i.op(1));
 207         const varidx & i2 = ex_to_varidx(i.op(2));
 208
 209         // A metric tensor with one covariant and one contravariant index gets
 210         // replaced by a delta tensor
 211         if (i1.is_covariant() != i2.is_covariant())
 212                 return delta_tensor(i1, i2);
 213
 214         // No further simplifications
 215         return i.hold();
 216 }
 217
 218 /** Automatic symbolic evaluation of an indexed Lorentz metric tensor. */
 219 ex minkmetric::eval_indexed(const basic & i) const
 220 {
 221         GINAC_ASSERT(is_of_type(i, indexed));
 222         GINAC_ASSERT(i.nops() == 3);
 223         GINAC_ASSERT(is_ex_of_type(i.op(0), minkmetric));
 224         GINAC_ASSERT(is_ex_of_type(i.op(1), varidx));
 225         GINAC_ASSERT(is_ex_of_type(i.op(2), varidx));
 226
 227         const varidx & i1 = ex_to_varidx(i.op(1));
 228         const varidx & i2 = ex_to_varidx(i.op(2));
 229
 230         // Numeric evaluation
 231         if (static_cast<const indexed &>(i).all_index_values_are(info_flags::nonnegint)) {
 232                 int n1 = ex_to_numeric(i1.get_value()).to_int(), n2 = ex_to_numeric(i2.get_value()).to_int();
 233                 if (n1 != n2)
 234                         return _ex0();
 235                 else if (n1 == 0)
 236                         return pos_sig ? _ex_1() : _ex1();
 237                 else
 238                         return pos_sig ? _ex1() : _ex_1();
 239         }
 240
 241         // Perform the usual evaluations of a metric tensor
 242         return inherited::eval_indexed(i);
 243 }
 244
 245 /** Automatic symbolic evaluation of an indexed epsilon tensor. */
 246 ex tensepsilon::eval_indexed(const basic & i) const
 247 {
 248         GINAC_ASSERT(is_of_type(i, indexed));
 249         GINAC_ASSERT(i.nops() > 1);
 250         GINAC_ASSERT(is_ex_of_type(i.op(0), tensepsilon));
 251
 252         // Convolutions are zero
 253         if (static_cast<const indexed &>(i).get_dummy_indices().size() != 0)
 254                 return _ex0();
 255
 256         // Numeric evaluation
 257         if (static_cast<const indexed &>(i).all_index_values_are(info_flags::nonnegint)) {
 258
 259                 // Get sign of index permutation (the indices should already be in
 260                 // a canonic order but we can't assume what exactly that order is)
 261                 std::vector<int> v;
 262                 v.reserve(i.nops() - 1);
 263                 for (unsigned j=1; j<i.nops(); j++)
 264                         v.push_back(ex_to_numeric(ex_to_idx(i.op(j)).get_value()).to_int());
 265                 int sign = permutation_sign(v);
 266
 267                 // In a Minkowski space, check for covariant indices
 268                 if (minkowski) {
 269                         for (unsigned j=1; j<i.nops(); j++) {
 270                                 const ex & x = i.op(j);
 271                                 if (!is_ex_of_type(x, varidx))
 272                                         throw(std::runtime_error("indices of epsilon tensor in Minkowski space must be of type varidx"));
 273                                 if (ex_to_varidx(x).is_covariant())
 274                                         if (ex_to_idx(x).get_value().is_zero())
 275                                                 sign = (pos_sig ? -sign : sign);
 276                                         else
 277                                                 sign = (pos_sig ? sign : -sign);
 278                         }
 279                 }
 280
 281                 return sign;
 282         }
 283
 284         // No further simplifications
 285         return i.hold();
 286 }
 287
 288 /** Contraction of an indexed delta tensor with something else. */
 289 bool tensdelta::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
 290 {
 291         GINAC_ASSERT(is_ex_of_type(*self, indexed));
 292         GINAC_ASSERT(is_ex_of_type(*other, indexed));
 293         GINAC_ASSERT(self->nops() == 3);
 294         GINAC_ASSERT(is_ex_of_type(self->op(0), tensdelta));
 295
 296         // Try to contract first index
 297         const idx *self_idx = &ex_to_idx(self->op(1));
 298         const idx *free_idx = &ex_to_idx(self->op(2));
 299         bool first_index_tried = false;
 300
 301 again:
 302         if (self_idx->is_symbolic()) {
 303                 for (int i=1; i<other->nops(); i++) {
 304                         const idx &other_idx = ex_to_idx(other->op(i));
 305                         if (is_dummy_pair(*self_idx, other_idx)) {
 306
 307                                 // Contraction found, remove delta tensor and substitute
 308                                 // index in second object
 309                                 *self = _ex1();
 310                                 *other = other->subs(other_idx == *free_idx);
 311                                 return true;
 312                         }
 313                 }
 314         }
 315
 316         if (!first_index_tried) {
 317
 318                 // No contraction with first index found, try second index
 319                 self_idx = &ex_to_idx(self->op(2));
 320                 free_idx = &ex_to_idx(self->op(1));
 321                 first_index_tried = true;
 322                 goto again;
 323         }
 324
 325         return false;
 326 }
 327
 328 /** Contraction of an indexed metric tensor with something else. */
 329 bool tensmetric::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
 330 {
 331         GINAC_ASSERT(is_ex_of_type(*self, indexed));
 332         GINAC_ASSERT(is_ex_of_type(*other, indexed));
 333         GINAC_ASSERT(self->nops() == 3);
 334         GINAC_ASSERT(is_ex_of_type(self->op(0), tensmetric));
 335
 336         // If contracting with the delta tensor, let the delta do it
 337         // (don't raise/lower delta indices)
 338         if (is_ex_of_type(other->op(0), tensdelta))
 339                 return false;
 340
 341         // Try to contract first index
 342         const idx *self_idx = &ex_to_idx(self->op(1));
 343         const idx *free_idx = &ex_to_idx(self->op(2));
 344         bool first_index_tried = false;
 345
 346 again:
 347         if (self_idx->is_symbolic()) {
 348                 for (int i=1; i<other->nops(); i++) {
 349                         const idx &other_idx = ex_to_idx(other->op(i));
 350                         if (is_dummy_pair(*self_idx, other_idx)) {
 351
 352                                 // Contraction found, remove metric tensor and substitute
 353                                 // index in second object
 354                                 *self = _ex1();
 355                                 *other = other->subs(other_idx == *free_idx);
 356                                 return true;
 357                         }
 358                 }
 359         }
 360
 361         if (!first_index_tried) {
 362
 363                 // No contraction with first index found, try second index
 364                 self_idx = &ex_to_idx(self->op(2));
 365                 free_idx = &ex_to_idx(self->op(1));
 366                 first_index_tried = true;
 367                 goto again;
 368         }
 369
 370         return false;
 371 }
 372
 373 //////////
 374 // global functions
 375 //////////
 376
 377 ex delta_tensor(const ex & i1, const ex & i2)
 378 {
 379         if (!is_ex_of_type(i1, idx) || !is_ex_of_type(i2, idx))
 380                 throw(std::invalid_argument("indices of delta tensor must be of type idx"));
 381
 382         return indexed(tensdelta(), indexed::symmetric, i1, i2);
 383 }
 384
 385 ex metric_tensor(const ex & i1, const ex & i2)
 386 {
 387         if (!is_ex_of_type(i1, varidx) || !is_ex_of_type(i2, varidx))
 388                 throw(std::invalid_argument("indices of metric tensor must be of type varidx"));
 389
 390         return indexed(tensmetric(), indexed::symmetric, i1, i2);
 391 }
 392
 393 ex lorentz_g(const ex & i1, const ex & i2, bool pos_sig)
 394 {
 395         if (!is_ex_of_type(i1, varidx) || !is_ex_of_type(i2, varidx))
 396                 throw(std::invalid_argument("indices of metric tensor must be of type varidx"));
 397
 398         return indexed(minkmetric(pos_sig), indexed::symmetric, i1, i2);
 399 }
 400
 401 ex epsilon_tensor(const ex & i1, const ex & i2)
 402 {
 403         if (!is_ex_of_type(i1, idx) || !is_ex_of_type(i2, idx))
 404                 throw(std::invalid_argument("indices of epsilon tensor must be of type idx"));
 405
 406         ex dim = ex_to_idx(i1).get_dim();
 407         if (!dim.is_equal(ex_to_idx(i2).get_dim()))
 408                 throw(std::invalid_argument("all indices of epsilon tensor must have the same dimension"));
 409         if (!ex_to_idx(i1).get_dim().is_equal(_ex2()))
 410                 throw(std::runtime_error("index dimension of epsilon tensor must match number of indices"));
 411
 412         return indexed(tensepsilon(), indexed::antisymmetric, i1, i2);
 413 }
 414
 415 ex epsilon_tensor(const ex & i1, const ex & i2, const ex & i3)
 416 {
 417         if (!is_ex_of_type(i1, idx) || !is_ex_of_type(i2, idx) || !is_ex_of_type(i3, idx))
 418                 throw(std::invalid_argument("indices of epsilon tensor must be of type idx"));
 419
 420         ex dim = ex_to_idx(i1).get_dim();
 421         if (!dim.is_equal(ex_to_idx(i2).get_dim()) || !dim.is_equal(ex_to_idx(i3).get_dim()))
 422                 throw(std::invalid_argument("all indices of epsilon tensor must have the same dimension"));
 423         if (!ex_to_idx(i1).get_dim().is_equal(_ex3()))
 424                 throw(std::runtime_error("index dimension of epsilon tensor must match number of indices"));
 425
 426         return indexed(tensepsilon(), indexed::antisymmetric, i1, i2, i3);
 427 }
 428
 429 ex lorentz_eps(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool pos_sig)
 430 {
 431         if (!is_ex_of_type(i1, varidx) || !is_ex_of_type(i2, varidx) || !is_ex_of_type(i3, varidx) || !is_ex_of_type(i4, varidx))
 432                 throw(std::invalid_argument("indices of Lorentz epsilon tensor must be of type varidx"));
 433
 434         ex dim = ex_to_idx(i1).get_dim();
 435         if (!dim.is_equal(ex_to_idx(i2).get_dim()) || !dim.is_equal(ex_to_idx(i3).get_dim()) || !dim.is_equal(ex_to_idx(i4).get_dim()))
 436                 throw(std::invalid_argument("all indices of epsilon tensor must have the same dimension"));
 437         if (!ex_to_idx(i1).get_dim().is_equal(_ex4()))
 438                 throw(std::runtime_error("index dimension of epsilon tensor must match number of indices"));
 439
 440         return indexed(tensepsilon(true, pos_sig), indexed::antisymmetric, i1, i2, i3, i4);
 441 }
 442
 443 } // namespace GiNaC