X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Ftensor.cpp;h=c215237d85f9d1225913b7fff4f87e59e0bc8134;hp=ccae243a6481377d36d460a147afed668a47231f;hb=ece407276fc0f844b2cfcd51838e1cbe83a2d70c;hpb=2e3eaaf5665eba012287220683995817cb371a13 diff --git a/ginac/tensor.cpp b/ginac/tensor.cpp index ccae243a..c215237d 100644 --- a/ginac/tensor.cpp +++ b/ginac/tensor.cpp @@ -3,7 +3,7 @@ * Implementation of GiNaC's special tensors. */ /* - * GiNaC Copyright (C) 1999-2002 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2006 Johannes Gutenberg University Mainz, Germany * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -17,7 +17,7 @@ * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ #include @@ -29,70 +29,72 @@ #include "indexed.h" #include "symmetry.h" #include "relational.h" +#include "operators.h" #include "lst.h" #include "numeric.h" #include "matrix.h" -#include "print.h" #include "archive.h" #include "utils.h" namespace GiNaC { 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) + +GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(tensdelta, tensor, + print_func(&tensdelta::do_print). + print_func(&tensdelta::do_print_latex)) + +GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(tensmetric, tensor, + print_func(&tensmetric::do_print). + print_func(&tensmetric::do_print)) + +GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(minkmetric, tensmetric, + print_func(&minkmetric::do_print). + print_func(&minkmetric::do_print_latex)) + +GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(spinmetric, tensmetric, + print_func(&spinmetric::do_print). + print_func(&spinmetric::do_print_latex)) + +GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(tensepsilon, tensor, + print_func(&tensepsilon::do_print). + print_func(&tensepsilon::do_print_latex)) ////////// -// default ctor, dtor, copy ctor, assignment operator and helpers +// constructors ////////// -DEFAULT_CTORS(tensor) -DEFAULT_CTORS(tensdelta) -DEFAULT_CTORS(tensmetric) -DEFAULT_COPY(spinmetric) -DEFAULT_DESTROY(spinmetric) -DEFAULT_DESTROY(minkmetric) -DEFAULT_DESTROY(tensepsilon) +tensor::tensor() : inherited(&tensor::tinfo_static) +{ + setflag(status_flags::evaluated | status_flags::expanded); +} + +DEFAULT_CTOR(tensdelta) +DEFAULT_CTOR(tensmetric) minkmetric::minkmetric() : pos_sig(false) { - tinfo_key = TINFO_minkmetric; + tinfo_key = &minkmetric::tinfo_static; } spinmetric::spinmetric() { - tinfo_key = TINFO_spinmetric; + tinfo_key = &spinmetric::tinfo_static; } minkmetric::minkmetric(bool ps) : pos_sig(ps) { - tinfo_key = TINFO_minkmetric; -} - -void minkmetric::copy(const minkmetric & other) -{ - inherited::copy(other); - pos_sig = other.pos_sig; + tinfo_key = &minkmetric::tinfo_static; } tensepsilon::tensepsilon() : minkowski(false), pos_sig(false) { - tinfo_key = TINFO_tensepsilon; + tinfo_key = &tensepsilon::tinfo_static; } tensepsilon::tensepsilon(bool mink, bool ps) : minkowski(mink), pos_sig(ps) { - tinfo_key = TINFO_tensepsilon; -} - -void tensepsilon::copy(const tensepsilon & other) -{ - inherited::copy(other); - minkowski = other.minkowski; - pos_sig = other.pos_sig; + tinfo_key = &tensepsilon::tinfo_static; } ////////// @@ -106,7 +108,7 @@ DEFAULT_ARCHIVING(spinmetric) DEFAULT_UNARCHIVE(minkmetric) DEFAULT_UNARCHIVE(tensepsilon) -minkmetric::minkmetric(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst) +minkmetric::minkmetric(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst) { n.find_bool("pos_sig", pos_sig); } @@ -117,7 +119,7 @@ void minkmetric::archive(archive_node &n) const n.add_bool("pos_sig", pos_sig); } -tensepsilon::tensepsilon(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst) +tensepsilon::tensepsilon(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst) { n.find_bool("minkowski", minkowski); n.find_bool("pos_sig", pos_sig); @@ -139,6 +141,22 @@ DEFAULT_COMPARE(tensdelta) DEFAULT_COMPARE(tensmetric) DEFAULT_COMPARE(spinmetric) +bool tensdelta::info(unsigned inf) const +{ + if(inf == info_flags::real) + return true; + + return false; +} + +bool tensmetric::info(unsigned inf) const +{ + if(inf == info_flags::real) + return true; + + return false; +} + int minkmetric::compare_same_type(const basic & other) const { GINAC_ASSERT(is_a(other)); @@ -150,6 +168,14 @@ int minkmetric::compare_same_type(const basic & other) const return inherited::compare_same_type(other); } +bool minkmetric::info(unsigned inf) const +{ + if(inf == info_flags::real) + return true; + + return false; +} + int tensepsilon::compare_same_type(const basic & other) const { GINAC_ASSERT(is_a(other)); @@ -163,6 +189,22 @@ int tensepsilon::compare_same_type(const basic & other) const return inherited::compare_same_type(other); } +bool tensepsilon::info(unsigned inf) const +{ + if(inf == info_flags::real) + return true; + + return false; +} + +bool spinmetric::info(unsigned inf) const +{ + if(inf == info_flags::real) + return true; + + return false; +} + DEFAULT_PRINT_LATEX(tensdelta, "delta", "\\delta") DEFAULT_PRINT(tensmetric, "g") DEFAULT_PRINT_LATEX(minkmetric, "eta", "\\eta") @@ -179,9 +221,24 @@ ex tensdelta::eval_indexed(const basic & i) const const idx & i1 = ex_to(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)) - return i1.get_dim(); + // The dimension of the indices must be equal, otherwise we use the minimal + // dimension + if (!i1.get_dim().is_equal(i2.get_dim())) { + ex min_dim = i1.minimal_dim(i2); + exmap m; + m[i1] = i1.replace_dim(min_dim); + m[i2] = i2.replace_dim(min_dim); + return i.subs(m, subs_options::no_pattern); + } + + // Trace of delta tensor is the (effective) dimension of the space + if (is_dummy_pair(i1, i2)) { + try { + return i1.minimal_dim(i2); + } catch (std::exception &e) { + return i.hold(); + } + } // Numeric evaluation if (static_cast(i).all_index_values_are(info_flags::integer)) { @@ -208,6 +265,16 @@ ex tensmetric::eval_indexed(const basic & i) const const varidx & i1 = ex_to(i.op(1)); const varidx & i2 = ex_to(i.op(2)); + // The dimension of the indices must be equal, otherwise we use the minimal + // dimension + if (!i1.get_dim().is_equal(i2.get_dim())) { + ex min_dim = i1.minimal_dim(i2); + exmap m; + m[i1] = i1.replace_dim(min_dim); + m[i2] = i2.replace_dim(min_dim); + return i.subs(m, subs_options::no_pattern); + } + // A metric tensor with one covariant and one contravariant index gets // replaced by a delta tensor if (i1.is_covariant() != i2.is_covariant()) @@ -293,15 +360,15 @@ ex tensepsilon::eval_indexed(const basic & i) const // a canonic order but we can't assume what exactly that order is) 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) { - for (unsigned j=1; j(x)) throw(std::runtime_error("indices of epsilon tensor in Minkowski space must be of type varidx")); if (ex_to(x).is_covariant()) if (ex_to(x).get_value().is_zero()) @@ -318,37 +385,37 @@ ex tensepsilon::eval_indexed(const basic & i) const return i.hold(); } -/** Contraction of an indexed delta tensor with something else. */ -bool tensdelta::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const +bool tensor::replace_contr_index(exvector::iterator self, exvector::iterator other) const { - GINAC_ASSERT(is_a(*self)); - GINAC_ASSERT(is_a(*other)); - GINAC_ASSERT(self->nops() == 3); - GINAC_ASSERT(is_a(self->op(0))); - - // Try to contract first index + // Try to contract the 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; again: if (self_idx->is_symbolic()) { - for (unsigned i=1; inops(); i++) { + for (size_t 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 - // index in second object - *self = _ex1; - *other = other->subs(other_idx == *free_idx); - return true; + // Contraction found, remove this tensor and substitute the + // index in the second object + try { + // minimal_dim() throws an exception when index dimensions are not comparable + ex min_dim = self_idx->minimal_dim(other_idx); + *other = other->subs(other_idx == free_idx->replace_dim(min_dim)); + *self = _ex1; // *other is assigned first because assigning *self invalidates free_idx + return true; + } catch (std::exception &e) { + return false; + } } } } if (!first_index_tried) { - // No contraction with first index found, try second index + // No contraction with the first index found, try the second index self_idx = &ex_to(self->op(2)); free_idx = &ex_to(self->op(1)); first_index_tried = true; @@ -358,6 +425,19 @@ again: return false; } +/** Contraction of an indexed delta tensor with something else. */ +bool tensdelta::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const +{ + GINAC_ASSERT(is_a(*self)); + GINAC_ASSERT(is_a(*other)); + GINAC_ASSERT(self->nops() == 3); + GINAC_ASSERT(is_a(self->op(0))); + + // Replace the dummy index with this tensor's other index and remove + // the tensor (this is valid for contractions with all other tensors) + return replace_contr_index(self, other); +} + /** Contraction of an indexed metric tensor with something else. */ bool tensmetric::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const { @@ -368,39 +448,12 @@ bool tensmetric::contract_with(exvector::iterator self, exvector::iterator other // 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)) + if (is_a(other->op(0))) 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; - -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 = _ex1; - *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; - goto again; - } - - return false; + // Replace the dummy index with this tensor's other index and remove + // the tensor + return replace_contr_index(self, other); } /** Contraction of an indexed spinor metric with something else. */ @@ -412,7 +465,7 @@ bool spinmetric::contract_with(exvector::iterator self, exvector::iterator other GINAC_ASSERT(is_a(self->op(0))); // Contractions between spinor metrics - if (is_ex_of_type(other->op(0), spinmetric)) { + if (is_a(other->op(0))) { 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)); @@ -445,7 +498,7 @@ bool spinmetric::contract_with(exvector::iterator self, exvector::iterator other // 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)) + if (is_a(other->op(0))) return false; // Try to contract first index @@ -456,14 +509,15 @@ bool spinmetric::contract_with(exvector::iterator self, exvector::iterator other again: if (self_idx->is_symbolic()) { - for (unsigned i=1; inops(); i++) { + for (size_t 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); + // index in second object (assign *self last because this + // invalidates free_idx) *other = other->subs(other_idx == *free_idx); + *self = (static_cast(self_idx)->is_covariant() ? sign : -sign); return true; } } @@ -488,67 +542,27 @@ bool tensepsilon::contract_with(exvector::iterator self, exvector::iterator othe GINAC_ASSERT(is_a(*self)); GINAC_ASSERT(is_a(*other)); GINAC_ASSERT(is_a(self->op(0))); - unsigned num = self->nops() - 1; + size_t num = self->nops() - 1; - if (is_ex_exactly_of_type(other->op(0), tensepsilon) && num+1 == other->nops()) { + if (is_exactly_a(other->op(0)) && num+1 == other->nops()) { // Contraction of two epsilon tensors is a determinant - ex dim = ex_to(self->op(1)).get_dim(); + bool variance = is_a(self->op(1)); matrix M(num, num); - for (int i=0; iop(i+1), other->op(j+1), pos_sig); - else + else if (variance) M(i, j) = metric_tensor(self->op(i+1), other->op(j+1)); + else + M(i, j) = delta_tensor(self->op(i+1), other->op(j+1)); } } int sign = minkowski ? -1 : 1; *self = sign * M.determinant().simplify_indexed(); *other = _ex1; return true; - - } else if (other->return_type() == return_types::commutative) { - -#if 0 - // This handles eps.i.j.k * p.j * p.k = 0 and related cases. - // Actually, simplify_indexed() can handle most of them on its own - // but one specific case that is not covered there is - // eps~mu.nu~i~j * p.mu * p~nu - // because of the difference in variance in the dummy indices mu - // and nu. Eventually, simplify_indexed() should be extended to - // handle this case, and this hack removed. - exvector c; - - // Handle all indices of the epsilon tensor - for (int i=0; iop(i+1); - - // Look whether there's a contraction with this index - exvector::const_iterator ait, aitend = v.end(); - for (ait = v.begin(); ait != aitend; ait++) { - if (ait == self) - continue; - if (is_a(*ait) && ait->return_type() == return_types::commutative && ex_to(*ait).has_dummy_index_for(idx) && ait->nops() == 2) { - - // Yes, did we already have another contraction with the same base expression? - ex base = ait->op(0); - if (std::find_if(c.begin(), c.end(), bind2nd(ex_is_equal(), base)) == c.end()) { - - // No, add the base expression to the list - c.push_back(base); - - } else { - - // Yes, the contraction is zero - *self = _ex0; - *other = _ex0; - return true; - } - } - } - } -#endif } return false; @@ -560,47 +574,52 @@ bool tensepsilon::contract_with(exvector::iterator self, exvector::iterator othe ex delta_tensor(const ex & i1, const ex & i2) { - if (!is_ex_of_type(i1, idx) || !is_ex_of_type(i2, idx)) + static ex delta = (new tensdelta)->setflag(status_flags::dynallocated); + + if (!is_a(i1) || !is_a(i2)) throw(std::invalid_argument("indices of delta tensor must be of type idx")); - return indexed(tensdelta(), sy_symm(), i1, i2); + return indexed(delta, symmetric2(), i1, i2); } ex metric_tensor(const ex & i1, const ex & i2) { - if (!is_ex_of_type(i1, varidx) || !is_ex_of_type(i2, varidx)) + static ex metric = (new tensmetric)->setflag(status_flags::dynallocated); + + if (!is_a(i1) || !is_a(i2)) throw(std::invalid_argument("indices of metric tensor must be of type varidx")); - ex dim = ex_to(i1).get_dim(); - if (!dim.is_equal(ex_to(i2).get_dim())) - throw(std::invalid_argument("all indices of metric tensor must have the same dimension")); - return indexed(tensmetric(), sy_symm(), i1, i2); + return indexed(metric, symmetric2(), i1, i2); } 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)) + static ex metric_neg = (new minkmetric(false))->setflag(status_flags::dynallocated); + static ex metric_pos = (new minkmetric(true))->setflag(status_flags::dynallocated); + + if (!is_a(i1) || !is_a(i2)) throw(std::invalid_argument("indices of metric tensor must be of type varidx")); - ex dim = ex_to(i1).get_dim(); - if (!dim.is_equal(ex_to(i2).get_dim())) - throw(std::invalid_argument("all indices of metric tensor must have the same dimension")); - return indexed(minkmetric(pos_sig), sy_symm(), i1, i2); + return indexed(pos_sig ? metric_pos : metric_neg, symmetric2(), i1, i2); } ex spinor_metric(const ex & i1, const ex & i2) { - if (!is_ex_of_type(i1, spinidx) || !is_ex_of_type(i2, spinidx)) + static ex metric = (new spinmetric)->setflag(status_flags::dynallocated); + + if (!is_a(i1) || !is_a(i2)) 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); + return indexed(metric, antisymmetric2(), i1, i2); } ex epsilon_tensor(const ex & i1, const ex & i2) { - if (!is_ex_of_type(i1, idx) || !is_ex_of_type(i2, idx)) + static ex epsilon = (new tensepsilon)->setflag(status_flags::dynallocated); + + if (!is_a(i1) || !is_a(i2)) throw(std::invalid_argument("indices of epsilon tensor must be of type idx")); ex dim = ex_to(i1).get_dim(); @@ -609,12 +628,17 @@ ex epsilon_tensor(const ex & i1, const ex & i2) 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(), sy_anti(), i1, i2); + if(is_a(i1.op(0))||is_a(i2.op(0))) + return indexed(epsilon, antisymmetric2(), i1, i2).hold(); + + return indexed(epsilon, antisymmetric2(), i1, i2); } 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)) + static ex epsilon = (new tensepsilon)->setflag(status_flags::dynallocated); + + if (!is_a(i1) || !is_a(i2) || !is_a(i3)) throw(std::invalid_argument("indices of epsilon tensor must be of type idx")); ex dim = ex_to(i1).get_dim(); @@ -623,12 +647,18 @@ ex epsilon_tensor(const ex & i1, const ex & i2, const ex & i3) 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(), sy_anti(), i1, i2, i3); + if(is_a(i1.op(0))||is_a(i2.op(0))||is_a(i3.op(0))) + return indexed(epsilon, antisymmetric3(), i1, i2, i3).hold(); + + return indexed(epsilon, antisymmetric3(), i1, i2, i3); } ex lorentz_eps(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)) + static ex epsilon_neg = (new tensepsilon(true, false))->setflag(status_flags::dynallocated); + static ex epsilon_pos = (new tensepsilon(true, true))->setflag(status_flags::dynallocated); + + if (!is_a(i1) || !is_a(i2) || !is_a(i3) || !is_a(i4)) throw(std::invalid_argument("indices of Lorentz epsilon tensor must be of type varidx")); ex dim = ex_to(i1).get_dim(); @@ -637,19 +667,10 @@ ex lorentz_eps(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool 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), sy_anti(), i1, i2, i3, i4); -} + if(is_a(i1.op(0))||is_a(i2.op(0))||is_a(i3.op(0))||is_a(i4.op(0))) + return indexed(pos_sig ? epsilon_pos : epsilon_neg, antisymmetric4(), i1, i2, i3, i4).hold(); -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); + return indexed(pos_sig ? epsilon_pos : epsilon_neg, antisymmetric4(), i1, i2, i3, i4); } } // namespace GiNaC