X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Ftensor.cpp;h=28a7d7bf61b7d8831ce8c3c312b4027668c49b73;hp=b53d6f7d0edb519a8ddee3ed6f5bd02855547971;hb=c4c7b40ea71dd20ecc8942f9cab50c4bb9766624;hpb=aef64ef1f63f7b4cb9236c5c189873e48c6c999b;ds=sidebyside diff --git a/ginac/tensor.cpp b/ginac/tensor.cpp index b53d6f7d..28a7d7bf 100644 --- a/ginac/tensor.cpp +++ b/ginac/tensor.cpp @@ -26,9 +26,12 @@ #include "tensor.h" #include "idx.h" #include "indexed.h" +#include "symmetry.h" #include "relational.h" #include "lst.h" #include "numeric.h" +#include "matrix.h" +#include "print.h" #include "archive.h" #include "utils.h" #include "debugmsg.h" @@ -39,6 +42,7 @@ GINAC_IMPLEMENT_REGISTERED_CLASS(tensor, basic) GINAC_IMPLEMENT_REGISTERED_CLASS(tensdelta, tensor) GINAC_IMPLEMENT_REGISTERED_CLASS(tensmetric, tensor) GINAC_IMPLEMENT_REGISTERED_CLASS(minkmetric, tensmetric) +GINAC_IMPLEMENT_REGISTERED_CLASS(spinmetric, tensmetric) GINAC_IMPLEMENT_REGISTERED_CLASS(tensepsilon, tensor) ////////// @@ -53,6 +57,8 @@ tensor::tensor(unsigned ti) : inherited(ti) DEFAULT_CTORS(tensor) DEFAULT_CTORS(tensdelta) DEFAULT_CTORS(tensmetric) +DEFAULT_COPY(spinmetric) +DEFAULT_DESTROY(spinmetric) DEFAULT_DESTROY(minkmetric) DEFAULT_DESTROY(tensepsilon) @@ -62,6 +68,12 @@ minkmetric::minkmetric() : pos_sig(false) tinfo_key = TINFO_minkmetric; } +spinmetric::spinmetric() +{ + debugmsg("spinmetric default constructor", LOGLEVEL_CONSTRUCT); + tinfo_key = TINFO_spinmetric; +} + minkmetric::minkmetric(bool ps) : pos_sig(ps) { debugmsg("minkmetric constructor from bool", LOGLEVEL_CONSTRUCT); @@ -100,6 +112,7 @@ void tensepsilon::copy(const tensepsilon & other) DEFAULT_ARCHIVING(tensor) DEFAULT_ARCHIVING(tensdelta) DEFAULT_ARCHIVING(tensmetric) +DEFAULT_ARCHIVING(spinmetric) DEFAULT_UNARCHIVE(minkmetric) DEFAULT_UNARCHIVE(tensepsilon) @@ -130,12 +143,13 @@ void tensepsilon::archive(archive_node &n) const } ////////// -// functions overriding virtual functions from bases classes +// functions overriding virtual functions from base classes ////////// DEFAULT_COMPARE(tensor) DEFAULT_COMPARE(tensdelta) DEFAULT_COMPARE(tensmetric) +DEFAULT_COMPARE(spinmetric) int minkmetric::compare_same_type(const basic & other) const { @@ -161,22 +175,11 @@ int tensepsilon::compare_same_type(const basic & other) const return inherited::compare_same_type(other); } -DEFAULT_PRINT(tensdelta, "delta") +DEFAULT_PRINT_LATEX(tensdelta, "delta", "\\delta") DEFAULT_PRINT(tensmetric, "g") -DEFAULT_PRINT(minkmetric, "eta") -DEFAULT_PRINT(tensepsilon, "eps") - -ex tensor::subs(const lst & ls, const lst & lr) const -{ - GINAC_ASSERT(ls.nops() == lr.nops()); - - for (unsigned i=0; i(i.op(1)); + const idx & i2 = ex_to(i.op(2)); // Trace of delta tensor is the dimension of the space if (is_dummy_pair(i1, i2)) @@ -194,7 +197,7 @@ ex tensdelta::eval_indexed(const basic & i) const // Numeric evaluation if (static_cast(i).all_index_values_are(info_flags::integer)) { - int n1 = ex_to_numeric(i1.get_value()).to_int(), n2 = ex_to_numeric(i2.get_value()).to_int(); + int n1 = ex_to(i1.get_value()).to_int(), n2 = ex_to(i2.get_value()).to_int(); if (n1 == n2) return _ex1(); else @@ -214,8 +217,8 @@ ex tensmetric::eval_indexed(const basic & i) const GINAC_ASSERT(is_ex_of_type(i.op(1), varidx)); GINAC_ASSERT(is_ex_of_type(i.op(2), varidx)); - const varidx & i1 = ex_to_varidx(i.op(1)); - const varidx & i2 = ex_to_varidx(i.op(2)); + const varidx & i1 = ex_to(i.op(1)); + const varidx & i2 = ex_to(i.op(2)); // A metric tensor with one covariant and one contravariant index gets // replaced by a delta tensor @@ -235,12 +238,12 @@ ex minkmetric::eval_indexed(const basic & i) const GINAC_ASSERT(is_ex_of_type(i.op(1), varidx)); GINAC_ASSERT(is_ex_of_type(i.op(2), varidx)); - const varidx & i1 = ex_to_varidx(i.op(1)); - const varidx & i2 = ex_to_varidx(i.op(2)); + const varidx & i1 = ex_to(i.op(1)); + const varidx & i2 = ex_to(i.op(2)); // Numeric evaluation if (static_cast(i).all_index_values_are(info_flags::nonnegint)) { - int n1 = ex_to_numeric(i1.get_value()).to_int(), n2 = ex_to_numeric(i2.get_value()).to_int(); + int n1 = ex_to(i1.get_value()).to_int(), n2 = ex_to(i2.get_value()).to_int(); if (n1 != n2) return _ex0(); else if (n1 == 0) @@ -253,6 +256,37 @@ ex minkmetric::eval_indexed(const basic & i) const return inherited::eval_indexed(i); } +/** Automatic symbolic evaluation of an indexed metric tensor. */ +ex spinmetric::eval_indexed(const basic & i) const +{ + GINAC_ASSERT(is_of_type(i, indexed)); + GINAC_ASSERT(i.nops() == 3); + GINAC_ASSERT(is_ex_of_type(i.op(0), spinmetric)); + GINAC_ASSERT(is_ex_of_type(i.op(1), spinidx)); + GINAC_ASSERT(is_ex_of_type(i.op(2), spinidx)); + + const spinidx & i1 = ex_to(i.op(1)); + const spinidx & i2 = ex_to(i.op(2)); + + // Convolutions are zero + if (!(static_cast(i).get_dummy_indices().empty())) + return _ex0(); + + // Numeric evaluation + if (static_cast(i).all_index_values_are(info_flags::nonnegint)) { + int n1 = ex_to(i1.get_value()).to_int(), n2 = ex_to(i2.get_value()).to_int(); + if (n1 == n2) + return _ex0(); + else if (n1 < n2) + return _ex1(); + else + return _ex_1(); + } + + // No further simplifications + return i.hold(); +} + /** Automatic symbolic evaluation of an indexed epsilon tensor. */ ex tensepsilon::eval_indexed(const basic & i) const { @@ -261,7 +295,7 @@ ex tensepsilon::eval_indexed(const basic & i) const GINAC_ASSERT(is_ex_of_type(i.op(0), tensepsilon)); // Convolutions are zero - if (static_cast(i).get_dummy_indices().size() != 0) + if (!(static_cast(i).get_dummy_indices().empty())) return _ex0(); // Numeric evaluation @@ -272,8 +306,8 @@ ex tensepsilon::eval_indexed(const basic & i) const std::vector v; v.reserve(i.nops() - 1); for (unsigned j=1; j(ex_to(i.op(j)).get_value()).to_int()); + int sign = permutation_sign(v.begin(), v.end()); // In a Minkowski space, check for covariant indices if (minkowski) { @@ -281,8 +315,8 @@ ex tensepsilon::eval_indexed(const basic & i) const const ex & x = i.op(j); if (!is_ex_of_type(x, varidx)) throw(std::runtime_error("indices of epsilon tensor in Minkowski space must be of type varidx")); - if (ex_to_varidx(x).is_covariant()) - if (ex_to_idx(x).get_value().is_zero()) + if (ex_to(x).is_covariant()) + if (ex_to(x).get_value().is_zero()) sign = (pos_sig ? -sign : sign); else sign = (pos_sig ? sign : -sign); @@ -305,14 +339,14 @@ bool tensdelta::contract_with(exvector::iterator self, exvector::iterator other, GINAC_ASSERT(is_ex_of_type(self->op(0), tensdelta)); // Try to contract first index - const idx *self_idx = &ex_to_idx(self->op(1)); - const idx *free_idx = &ex_to_idx(self->op(2)); + const idx *self_idx = &ex_to(self->op(1)); + const idx *free_idx = &ex_to(self->op(2)); bool first_index_tried = false; again: if (self_idx->is_symbolic()) { - for (int i=1; inops(); i++) { - const idx &other_idx = ex_to_idx(other->op(i)); + for (unsigned i=1; inops(); i++) { + const idx &other_idx = ex_to(other->op(i)); if (is_dummy_pair(*self_idx, other_idx)) { // Contraction found, remove delta tensor and substitute @@ -327,8 +361,8 @@ again: if (!first_index_tried) { // No contraction with first index found, try second index - self_idx = &ex_to_idx(self->op(2)); - free_idx = &ex_to_idx(self->op(1)); + self_idx = &ex_to(self->op(2)); + free_idx = &ex_to(self->op(1)); first_index_tried = true; goto again; } @@ -350,14 +384,14 @@ bool tensmetric::contract_with(exvector::iterator self, exvector::iterator other return false; // Try to contract first index - const idx *self_idx = &ex_to_idx(self->op(1)); - const idx *free_idx = &ex_to_idx(self->op(2)); + const idx *self_idx = &ex_to(self->op(1)); + const idx *free_idx = &ex_to(self->op(2)); bool first_index_tried = false; again: if (self_idx->is_symbolic()) { - for (int i=1; inops(); i++) { - const idx &other_idx = ex_to_idx(other->op(i)); + for (unsigned i=1; inops(); i++) { + const idx &other_idx = ex_to(other->op(i)); if (is_dummy_pair(*self_idx, other_idx)) { // Contraction found, remove metric tensor and substitute @@ -372,15 +406,119 @@ again: if (!first_index_tried) { // No contraction with first index found, try second index - self_idx = &ex_to_idx(self->op(2)); - free_idx = &ex_to_idx(self->op(1)); + self_idx = &ex_to(self->op(2)); + free_idx = &ex_to(self->op(1)); + first_index_tried = true; + goto again; + } + + return false; +} + +/** Contraction of an indexed spinor metric with something else. */ +bool spinmetric::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const +{ + GINAC_ASSERT(is_ex_of_type(*self, indexed)); + GINAC_ASSERT(is_ex_of_type(*other, indexed)); + GINAC_ASSERT(self->nops() == 3); + GINAC_ASSERT(is_ex_of_type(self->op(0), spinmetric)); + + // Contractions between spinor metrics + if (is_ex_of_type(other->op(0), spinmetric)) { + const idx &self_i1 = ex_to(self->op(1)); + const idx &self_i2 = ex_to(self->op(2)); + const idx &other_i1 = ex_to(other->op(1)); + const idx &other_i2 = ex_to(other->op(2)); + + if (is_dummy_pair(self_i1, other_i1)) { + if (is_dummy_pair(self_i2, other_i2)) + *self = _ex2(); + else + *self = delta_tensor(self_i2, other_i2); + *other = _ex1(); + return true; + } else if (is_dummy_pair(self_i1, other_i2)) { + if (is_dummy_pair(self_i2, other_i1)) + *self = _ex_2(); + else + *self = -delta_tensor(self_i2, other_i1); + *other = _ex1(); + return true; + } else if (is_dummy_pair(self_i2, other_i1)) { + *self = -delta_tensor(self_i1, other_i2); + *other = _ex1(); + return true; + } else if (is_dummy_pair(self_i2, other_i2)) { + *self = delta_tensor(self_i1, other_i1); + *other = _ex1(); + return true; + } + } + + // If contracting with the delta tensor, let the delta do it + // (don't raise/lower delta indices) + if (is_ex_of_type(other->op(0), tensdelta)) + return false; + + // Try to contract first index + const idx *self_idx = &ex_to(self->op(1)); + const idx *free_idx = &ex_to(self->op(2)); + bool first_index_tried = false; + int sign = 1; + +again: + if (self_idx->is_symbolic()) { + for (unsigned i=1; inops(); i++) { + const idx &other_idx = ex_to(other->op(i)); + if (is_dummy_pair(*self_idx, other_idx)) { + + // Contraction found, remove metric tensor and substitute + // index in second object + *self = (static_cast(self_idx)->is_covariant() ? sign : -sign); + *other = other->subs(other_idx == *free_idx); + return true; + } + } + } + + if (!first_index_tried) { + + // No contraction with first index found, try second index + self_idx = &ex_to(self->op(2)); + free_idx = &ex_to(self->op(1)); first_index_tried = true; + sign = -sign; goto again; } return false; } +/** Contraction of epsilon tensor with something else. */ +bool tensepsilon::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const +{ + GINAC_ASSERT(is_ex_of_type(*self, indexed)); + GINAC_ASSERT(is_ex_of_type(*other, indexed)); + GINAC_ASSERT(is_ex_of_type(self->op(0), spinmetric)); + unsigned num = self->nops() - 1; + + if (is_ex_exactly_of_type(other->op(0), tensepsilon) && num+1 == other->nops()) { + + // Contraction of two epsilon tensors is a determinant + ex dim = ex_to(self->op(1)).get_dim(); + matrix M(num, num); + for (int i=0; iop(i+1), other->op(j+1)); + int sign = minkowski ? -1 : 1; + *self = sign * M.determinant().simplify_indexed(); + *other = _ex1(); + return true; + } + + return false; +} + ////////// // global functions ////////// @@ -390,7 +528,7 @@ ex delta_tensor(const ex & i1, const ex & i2) if (!is_ex_of_type(i1, idx) || !is_ex_of_type(i2, idx)) throw(std::invalid_argument("indices of delta tensor must be of type idx")); - return indexed(tensdelta(), indexed::symmetric, i1, i2); + return indexed(tensdelta(), sy_symm(), i1, i2); } ex metric_tensor(const ex & i1, const ex & i2) @@ -398,7 +536,7 @@ ex metric_tensor(const ex & i1, const ex & i2) if (!is_ex_of_type(i1, varidx) || !is_ex_of_type(i2, varidx)) throw(std::invalid_argument("indices of metric tensor must be of type varidx")); - return indexed(tensmetric(), indexed::symmetric, i1, i2); + return indexed(tensmetric(), sy_symm(), i1, i2); } ex lorentz_g(const ex & i1, const ex & i2, bool pos_sig) @@ -406,7 +544,17 @@ ex lorentz_g(const ex & i1, const ex & i2, bool pos_sig) if (!is_ex_of_type(i1, varidx) || !is_ex_of_type(i2, varidx)) throw(std::invalid_argument("indices of metric tensor must be of type varidx")); - return indexed(minkmetric(pos_sig), indexed::symmetric, i1, i2); + return indexed(minkmetric(pos_sig), sy_symm(), i1, i2); +} + +ex spinor_metric(const ex & i1, const ex & i2) +{ + if (!is_ex_of_type(i1, spinidx) || !is_ex_of_type(i2, spinidx)) + throw(std::invalid_argument("indices of spinor metric must be of type spinidx")); + if (!ex_to(i1).get_dim().is_equal(2) || !ex_to(i2).get_dim().is_equal(2)) + throw(std::runtime_error("index dimension for spinor metric must be 2")); + + return indexed(spinmetric(), sy_anti(), i1, i2); } ex epsilon_tensor(const ex & i1, const ex & i2) @@ -414,13 +562,13 @@ ex epsilon_tensor(const ex & i1, const ex & i2) if (!is_ex_of_type(i1, idx) || !is_ex_of_type(i2, idx)) throw(std::invalid_argument("indices of epsilon tensor must be of type idx")); - ex dim = ex_to_idx(i1).get_dim(); - if (!dim.is_equal(ex_to_idx(i2).get_dim())) + ex dim = ex_to(i1).get_dim(); + if (!dim.is_equal(ex_to(i2).get_dim())) throw(std::invalid_argument("all indices of epsilon tensor must have the same dimension")); - if (!ex_to_idx(i1).get_dim().is_equal(_ex2())) + if (!ex_to(i1).get_dim().is_equal(_ex2())) throw(std::runtime_error("index dimension of epsilon tensor must match number of indices")); - return indexed(tensepsilon(), indexed::antisymmetric, i1, i2); + return indexed(tensepsilon(), sy_anti(), i1, i2); } ex epsilon_tensor(const ex & i1, const ex & i2, const ex & i3) @@ -428,13 +576,13 @@ ex epsilon_tensor(const ex & i1, const ex & i2, const ex & i3) if (!is_ex_of_type(i1, idx) || !is_ex_of_type(i2, idx) || !is_ex_of_type(i3, idx)) throw(std::invalid_argument("indices of epsilon tensor must be of type idx")); - ex dim = ex_to_idx(i1).get_dim(); - if (!dim.is_equal(ex_to_idx(i2).get_dim()) || !dim.is_equal(ex_to_idx(i3).get_dim())) + ex dim = ex_to(i1).get_dim(); + if (!dim.is_equal(ex_to(i2).get_dim()) || !dim.is_equal(ex_to(i3).get_dim())) throw(std::invalid_argument("all indices of epsilon tensor must have the same dimension")); - if (!ex_to_idx(i1).get_dim().is_equal(_ex3())) + if (!ex_to(i1).get_dim().is_equal(_ex3())) throw(std::runtime_error("index dimension of epsilon tensor must match number of indices")); - return indexed(tensepsilon(), indexed::antisymmetric, i1, i2, i3); + return indexed(tensepsilon(), sy_anti(), i1, i2, i3); } ex lorentz_eps(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool pos_sig) @@ -442,13 +590,25 @@ ex lorentz_eps(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool 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)) throw(std::invalid_argument("indices of Lorentz epsilon tensor must be of type varidx")); - ex dim = ex_to_idx(i1).get_dim(); - 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())) + ex dim = ex_to(i1).get_dim(); + if (!dim.is_equal(ex_to(i2).get_dim()) || !dim.is_equal(ex_to(i3).get_dim()) || !dim.is_equal(ex_to(i4).get_dim())) throw(std::invalid_argument("all indices of epsilon tensor must have the same dimension")); - if (!ex_to_idx(i1).get_dim().is_equal(_ex4())) + if (!ex_to(i1).get_dim().is_equal(_ex4())) throw(std::runtime_error("index dimension of epsilon tensor must match number of indices")); - return indexed(tensepsilon(true, pos_sig), indexed::antisymmetric, i1, i2, i3, i4); + return indexed(tensepsilon(true, pos_sig), sy_anti(), i1, i2, i3, i4); +} + +ex eps0123(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool pos_sig) +{ + 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)) + throw(std::invalid_argument("indices of epsilon tensor must be of type varidx")); + + ex dim = ex_to(i1).get_dim(); + if (dim.is_equal(4)) + return lorentz_eps(i1, i2, i3, i4, pos_sig); + else + return indexed(tensepsilon(true, pos_sig), sy_anti(), i1, i2, i3, i4); } } // namespace GiNaC