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
  25 #include "tensor.h"
  26 #include "idx.h"
  27 #include "indexed.h"
  28 #include "relational.h"
  29 #include "numeric.h"
  30 #include "archive.h"
  31 #include "utils.h"
  32 #include "debugmsg.h"
  33
  34 namespace GiNaC {
  35
  36 GINAC_IMPLEMENT_REGISTERED_CLASS(tensor, basic)
  37 GINAC_IMPLEMENT_REGISTERED_CLASS(tensdelta, tensor)
  38 GINAC_IMPLEMENT_REGISTERED_CLASS(tensmetric, tensor)
  39 GINAC_IMPLEMENT_REGISTERED_CLASS(minkmetric, tensmetric)
  40 GINAC_IMPLEMENT_REGISTERED_CLASS(tensepsilon, tensor)
  41
  42 //////////
  43 // default constructor, destructor, copy constructor assignment operator and helpers
  44 //////////
  45
  46 #define DEFAULT_CTORS(classname) \
  47 classname::classname() : inherited(TINFO_##classname) \
  48 { \
  49         debugmsg(#classname " default constructor", LOGLEVEL_CONSTRUCT); \
  50 } \
  51 void classname::copy(const classname & other) \
  52 { \
  53         inherited::copy(other); \
  54 } \
  55 void classname::destroy(bool call_parent) \
  56 { \
  57         if (call_parent) \
  58                 inherited::destroy(call_parent); \
  59 }
  60
  61 tensor::tensor(unsigned ti) : inherited(ti)
  62 {
  63         debugmsg("tensor constructor from unsigned", LOGLEVEL_CONSTRUCT); \
  64 }
  65
  66 DEFAULT_CTORS(tensor)
  67 DEFAULT_CTORS(tensdelta)
  68 DEFAULT_CTORS(tensmetric)
  69 DEFAULT_CTORS(tensepsilon)
  70
  71 minkmetric::minkmetric() : pos_sig(false)
  72 {
  73         debugmsg("minkmetric default constructor", LOGLEVEL_CONSTRUCT);
  74         tinfo_key = TINFO_minkmetric;
  75 }
  76
  77 minkmetric::minkmetric(bool ps) : pos_sig(ps)
  78 {
  79         debugmsg("minkmetric constructor from bool", LOGLEVEL_CONSTRUCT);
  80         tinfo_key = TINFO_minkmetric;
  81 }
  82
  83 void minkmetric::copy(const minkmetric & other)
  84 {
  85         inherited::copy(other);
  86         pos_sig = other.pos_sig;
  87 }
  88
  89 void minkmetric::destroy(bool call_parent)
  90 {
  91         if (call_parent)
  92                 inherited::destroy(call_parent);
  93 }
  94
  95 //////////
  96 // archiving
  97 //////////
  98
  99 #define DEFAULT_ARCHIVING(classname) \
 100 classname::classname(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst) \
 101 { \
 102         debugmsg(#classname " constructor from archive_node", LOGLEVEL_CONSTRUCT); \
 103 } \
 104 ex classname::unarchive(const archive_node &n, const lst &sym_lst) \
 105 { \
 106         return (new classname(n, sym_lst))->setflag(status_flags::dynallocated); \
 107 } \
 108 void classname::archive(archive_node &n) const \
 109 { \
 110         inherited::archive(n); \
 111 }
 112
 113 DEFAULT_ARCHIVING(tensor)
 114 DEFAULT_ARCHIVING(tensdelta)
 115 DEFAULT_ARCHIVING(tensmetric)
 116 DEFAULT_ARCHIVING(tensepsilon)
 117
 118 minkmetric::minkmetric(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
 119 {
 120         debugmsg("minkmetric constructor from archive_node", LOGLEVEL_CONSTRUCT);
 121         n.find_bool("pos_sig", pos_sig);
 122 }
 123
 124 ex minkmetric::unarchive(const archive_node &n, const lst &sym_lst)
 125 {
 126         return (new minkmetric(n, sym_lst))->setflag(status_flags::dynallocated);
 127 }
 128
 129 void minkmetric::archive(archive_node &n) const
 130 {
 131         inherited::archive(n);
 132         n.add_bool("pos_sig", pos_sig);
 133 }
 134
 135 //////////
 136 // functions overriding virtual functions from bases classes
 137 //////////
 138
 139 #define DEFAULT_COMPARE(classname) \
 140 int classname::compare_same_type(const basic & other) const \
 141 { \
 142         /* by default, two tensors of the same class are always identical */ \
 143         return 0; \
 144 }
 145
 146 DEFAULT_COMPARE(tensor)
 147 DEFAULT_COMPARE(tensdelta)
 148 DEFAULT_COMPARE(tensmetric)
 149 DEFAULT_COMPARE(tensepsilon)
 150
 151 int minkmetric::compare_same_type(const basic & other) const
 152 {
 153         GINAC_ASSERT(is_of_type(other, minkmetric));
 154         const minkmetric &o = static_cast<const minkmetric &>(other);
 155
 156         if (pos_sig != o.pos_sig)
 157                 return pos_sig ? -1 : 1;
 158         else
 159                 return inherited::compare_same_type(other);
 160 }
 161
 162 void tensdelta::print(std::ostream & os, unsigned upper_precedence) const
 163 {
 164         debugmsg("tensdelta print",LOGLEVEL_PRINT);
 165         os << "delta";
 166 }
 167
 168 void tensmetric::print(std::ostream & os, unsigned upper_precedence) const
 169 {
 170         debugmsg("tensmetric print",LOGLEVEL_PRINT);
 171         os << "g";
 172 }
 173
 174 void minkmetric::print(std::ostream & os, unsigned upper_precedence) const
 175 {
 176         debugmsg("minkmetric print",LOGLEVEL_PRINT);
 177         os << "eta";
 178 }
 179
 180 void tensepsilon::print(std::ostream & os, unsigned upper_precedence) const
 181 {
 182         debugmsg("tensepsilon print",LOGLEVEL_PRINT);
 183         os << "eps";
 184 }
 185
 186 /** Automatic symbolic evaluation of an indexed delta tensor. */
 187 ex tensdelta::eval_indexed(const basic & i) const
 188 {
 189         GINAC_ASSERT(is_of_type(i, indexed));
 190         GINAC_ASSERT(i.nops() == 3);
 191         GINAC_ASSERT(is_ex_of_type(i.op(0), tensdelta));
 192
 193         const idx & i1 = ex_to_idx(i.op(1));
 194         const idx & i2 = ex_to_idx(i.op(2));
 195
 196         // Trace of delta tensor is the dimension of the space
 197         if (is_dummy_pair(i1, i2))
 198                 return i1.get_dim();
 199
 200         // No further simplifications
 201         return i.hold();
 202 }
 203
 204 /** Automatic symbolic evaluation of an indexed metric tensor. */
 205 ex tensmetric::eval_indexed(const basic & i) const
 206 {
 207         GINAC_ASSERT(is_of_type(i, indexed));
 208         GINAC_ASSERT(i.nops() == 3);
 209         GINAC_ASSERT(is_ex_of_type(i.op(0), tensmetric));
 210         GINAC_ASSERT(is_ex_of_type(i.op(1), varidx));
 211         GINAC_ASSERT(is_ex_of_type(i.op(2), varidx));
 212
 213         const varidx & i1 = ex_to_varidx(i.op(1));
 214         const varidx & i2 = ex_to_varidx(i.op(2));
 215
 216         // A metric tensor with one covariant and one contravariant index gets
 217         // replaced by a delta tensor
 218         if (i1.is_covariant() != i2.is_covariant())
 219                 return delta_tensor(i1, i2);
 220
 221         // No further simplifications
 222         return i.hold();
 223 }
 224
 225 /** Automatic symbolic evaluation of an indexed Lorentz metric tensor. */
 226 ex minkmetric::eval_indexed(const basic & i) const
 227 {
 228         GINAC_ASSERT(is_of_type(i, indexed));
 229         GINAC_ASSERT(i.nops() == 3);
 230         GINAC_ASSERT(is_ex_of_type(i.op(0), minkmetric));
 231         GINAC_ASSERT(is_ex_of_type(i.op(1), varidx));
 232         GINAC_ASSERT(is_ex_of_type(i.op(2), varidx));
 233
 234         const varidx & i1 = ex_to_varidx(i.op(1));
 235         const varidx & i2 = ex_to_varidx(i.op(2));
 236
 237         // Numeric evaluation
 238         if (static_cast<const indexed &>(i).all_index_values_are(info_flags::nonnegint)) {
 239                 int n1 = ex_to_numeric(i1.get_value()).to_int(), n2 = ex_to_numeric(i2.get_value()).to_int();
 240                 if (n1 != n2)
 241                         return _ex0();
 242                 else if (n1 == 0)
 243                         return pos_sig ? _ex_1() : _ex1();
 244                 else
 245                         return pos_sig ? _ex1() : _ex_1();
 246         }
 247
 248         // Perform the usual evaluations of a metric tensor
 249         return inherited::eval_indexed(i);
 250 }
 251
 252 /** Contraction of an indexed delta tensor with something else. */
 253 bool tensdelta::contract_with(ex & self, ex & other) const
 254 {
 255         GINAC_ASSERT(is_ex_of_type(self, indexed));
 256         GINAC_ASSERT(is_ex_of_type(other, indexed));
 257         GINAC_ASSERT(self.nops() == 3);
 258         GINAC_ASSERT(is_ex_of_type(self.op(0), tensdelta));
 259
 260         // Try to contract first index
 261         const idx *self_idx = &ex_to_idx(self.op(1));
 262         const idx *free_idx = &ex_to_idx(self.op(2));
 263         bool first_index_tried = false;
 264
 265 again:
 266         if (self_idx->is_symbolic()) {
 267                 for (int i=1; i<other.nops(); i++) {
 268                         const idx &other_idx = ex_to_idx(other.op(i));
 269                         if (is_dummy_pair(*self_idx, other_idx)) {
 270
 271                                 // Contraction found, remove delta tensor and substitute
 272                                 // index in second object
 273                                 self = _ex1();
 274                                 other = other.subs(other_idx == *free_idx);
 275                                 return true;
 276                         }
 277                 }
 278         }
 279
 280         if (!first_index_tried) {
 281
 282                 // No contraction with first index found, try second index
 283                 self_idx = &ex_to_idx(self.op(2));
 284                 free_idx = &ex_to_idx(self.op(1));
 285                 first_index_tried = true;
 286                 goto again;
 287         }
 288
 289         return false;
 290 }
 291
 292 /** Contraction of an indexed metric tensor with something else. */
 293 bool tensmetric::contract_with(ex & self, ex & other) const
 294 {
 295         GINAC_ASSERT(is_ex_of_type(self, indexed));
 296         GINAC_ASSERT(is_ex_of_type(other, indexed));
 297         GINAC_ASSERT(self.nops() == 3);
 298         GINAC_ASSERT(is_ex_of_type(self.op(0), tensmetric));
 299
 300         // If contracting with the delta tensor, let the delta do it
 301         // (don't raise/lower delta indices)
 302         if (is_ex_exactly_of_type(other.op(0), tensdelta))
 303                 return false;
 304
 305         // Try to contract first index
 306         const idx *self_idx = &ex_to_idx(self.op(1));
 307         const idx *free_idx = &ex_to_idx(self.op(2));
 308         bool first_index_tried = false;
 309
 310 again:
 311         if (self_idx->is_symbolic()) {
 312                 for (int i=1; i<other.nops(); i++) {
 313                         const idx &other_idx = ex_to_idx(other.op(i));
 314                         if (is_dummy_pair(*self_idx, other_idx)) {
 315
 316                                 // Contraction found, remove metric tensor and substitute
 317                                 // index in second object
 318                                 self = _ex1();
 319                                 other = other.subs(other_idx == *free_idx);
 320                                 return true;
 321                         }
 322                 }
 323         }
 324
 325         if (!first_index_tried) {
 326
 327                 // No contraction with first index found, try second index
 328                 self_idx = &ex_to_idx(self.op(2));
 329                 free_idx = &ex_to_idx(self.op(1));
 330                 first_index_tried = true;
 331                 goto again;
 332         }
 333
 334         return false;
 335 }
 336
 337 //////////
 338 // global functions
 339 //////////
 340
 341 ex delta_tensor(const ex & i1, const ex & i2)
 342 {
 343         if (!is_ex_of_type(i1, idx) || !is_ex_of_type(i2, idx))
 344                 throw(std::invalid_argument("indices of delta tensor must be of type idx"));
 345
 346         return indexed(tensdelta(), indexed::symmetric, i1, i2);
 347 }
 348
 349 ex metric_tensor(const ex & i1, const ex & i2)
 350 {
 351         if (!is_ex_of_type(i1, varidx) || !is_ex_of_type(i2, varidx))
 352                 throw(std::invalid_argument("indices of metric tensor must be of type varidx"));
 353
 354         return indexed(tensmetric(), i1, i2);
 355 }
 356
 357 ex lorentz_g(const ex & i1, const ex & i2, bool pos_sig)
 358 {
 359         if (!is_ex_of_type(i1, varidx) || !is_ex_of_type(i2, varidx))
 360                 throw(std::invalid_argument("indices of metric tensor must be of type varidx"));
 361
 362         return indexed(minkmetric(pos_sig), indexed::symmetric, i1, i2);
 363 }
 364
 365 ex epsilon_tensor(const ex & i1, const ex & i2)
 366 {
 367         if (!is_ex_of_type(i1, idx) || !is_ex_of_type(i2, idx))
 368                 throw(std::invalid_argument("indices of epsilon tensor must be of type idx"));
 369         if (!ex_to_idx(i1).get_dim().is_equal(_ex2()) || !ex_to_idx(i2).get_dim().is_equal(_ex2()))
 370                 throw(std::invalid_argument("index dimension of epsilon tensor must match number of indices"));
 371
 372         return indexed(tensepsilon(), indexed::antisymmetric, i1, i2);
 373 }
 374
 375 } // namespace GiNaC