X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Ftensor.cpp;h=0c36a68002423d712f4d8afa5ec903716ffba7e7;hp=ccae243a6481377d36d460a147afed668a47231f;hb=b8b8cfbb72bfa59b01371f67f542914cf55f2ab9;hpb=2e3eaaf5665eba012287220683995817cb371a13 diff --git a/ginac/tensor.cpp b/ginac/tensor.cpp index ccae243a..0c36a680 100644 --- a/ginac/tensor.cpp +++ b/ginac/tensor.cpp @@ -179,9 +179,14 @@ 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(); + // 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 +213,13 @@ 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); + return i.subs(lst(i1 == i1.replace_dim(min_dim), i2 == i2.replace_dim(min_dim))); + } + // A metric tensor with one covariant and one contravariant index gets // replaced by a delta tensor if (i1.is_covariant() != i2.is_covariant()) @@ -318,15 +330,9 @@ 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; @@ -337,18 +343,24 @@ again: 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); + *self = _ex1; + *other = other->subs(other_idx == free_idx->replace_dim(min_dim)); + 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 +370,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 { @@ -371,36 +396,9 @@ bool tensmetric::contract_with(exvector::iterator self, exvector::iterator other 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; - -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 (this is valid for contractions with all other tensors) + return replace_contr_index(self, other); } /** Contraction of an indexed spinor metric with something else. */ @@ -507,48 +505,6 @@ bool tensepsilon::contract_with(exvector::iterator self, exvector::iterator othe *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; @@ -570,9 +526,6 @@ 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")); - 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); } @@ -581,9 +534,6 @@ 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")); - 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); } @@ -640,16 +590,4 @@ ex lorentz_eps(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool 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