[GiNaC-list] Covariant derivative(its long..)

[mdias at ift.unesp.br]

Hi,


I tried this very confuse code, and I am happy if somebody could help-me. I
tried to make covariant derivative derivated from a indexed. Use it like:
D(indexed(A,mu),nu)
 But the fatal errors in this code (among millions..):
i) simplify_indexed() crashs or doesnt work !!
II) get_free_indexes not work, so D^\alpha(A_beta)A^\beta A_\alpha
-D^\beta(A_alpha)A^\beta A_\alpha is not zero
 In the last post , Christian Bauer suggests "override is eval_indexed() (for
automatic evaluation) and
contract_with()/add_indexed()/scalar_mul_indexed() (for manual evaluation,
AKA simplify_indexed())" and showed the way: look color.cpp source.

Any suggestions are welcome

Marco

#include <iostream>
#include <ginac/ginac.h>


#define N 4
#define N1 8
using namespace std;
using namespace GiNaC;

const unsigned TINFO_covdiff = 0x42420001U;




class covdiff : public indexed {
GINAC_DECLARE_REGISTERED_CLASS(covdiff, indexed)
    public:
  covdiff(const ex  & b,const  ex & i1);
  covdiff(const ex & b, const ex & i1, const ex & i2);
  covdiff(const ex & b, const ex & i1, const ex & i2,const ex & i3);
  covdiff(const ex & b, const ex & i1, const ex & i2,const ex & i3, const ex & i4);
  covdiff(const ex & b, const lst & iv);
  size_t nops() const;
  ex op(size_t i) const; 
  ex &  let_op(size_t i);
  ex subs(const exmap & m, unsigned options) const;
  protected:
	void do_print(const print_context & c, unsigned level) const;
        ex symtree, arg,index; 
  lst seq;
};


GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(covdiff, indexed,
  print_func<print_context>(&covdiff::do_print))


  covdiff::covdiff(const ex & b, const ex & i1) : inherited(b,
i1),arg(b),index(i1), symtree(not_symmetric())
{
	tinfo_key = TINFO_covdiff;
}


covdiff::covdiff(const ex & b, const ex & i1, const ex & i2) : inherited(b, i1,
i2),arg(b), symtree(not_symmetric())
{
  seq =i1,i2;
        tinfo_key = TINFO_covdiff;
}

covdiff::covdiff(const ex & b, const ex & i1, const ex & i2, const ex & i3) :
inherited(b, i1, i2, i3),arg(b), symtree(not_symmetric())
{ seq=i1,i2,i3;
       tinfo_key = TINFO_covdiff;
}

covdiff::covdiff(const ex & b, const ex & i1, const ex & i2, const ex & i3,
const ex & i4) : inherited(b, i1, i2, i3, i4),arg(b), symtree(not_symmetric())
{seq = i1,i2,i3,i4;
        tinfo_key = TINFO_covdiff;
}
covdiff::covdiff(const ex & b, const lst & iv) : inherited(b),arg(b),
symtree(not_symmetric())
{
        seq=iv;
        tinfo_key = TINFO_covdiff;

}



//archiving
covdiff::covdiff(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst) 
{
  	if (!n.find_ex("symmetry", symtree, sym_lst)) {
		// GiNaC versions <= 0.9.0 had an unsigned "symmetry" property
		unsigned symm = 0;
		n.find_unsigned("symmetry", symm);
		switch (symm) {
			case 1:
				symtree = sy_symm();
				break;
			case 2:
				symtree = sy_anti();
				break;
			default:
				symtree = not_symmetric();
				break;
		}
		const_cast<symmetry &>(ex_to<symmetry>(symtree)).validate(seq.nops() - 1);
		}
}


void covdiff::archive(archive_node &n) const
{
	inherited::archive(n);
	n.add_ex("symmetry", symtree);
}


 
ex covdiff::unarchive(const archive_node &n, lst &sym_lst)
{
  return (new covdiff(n, sym_lst))->setflag(status_flags::dynallocated);
}

//comparation method


void covdiff::do_print(const print_context & c, unsigned level) const
{
if (is_a<print_tree>(c))
    covdiff::print(c,level);
if (seq.nops() > 1) {
  for( lst::const_iterator it=seq.begin();it != seq.end();++it)
    c.s<<"D"<<*it;
  c.s<<"("<<arg<<")";
}
 else 
  c.s << "D"<<index <<" ("<<arg<<")";
}



int covdiff::compare_same_type(const basic & other) const
{
GINAC_ASSERT(is_a<covdiff>(other));
	return inherited::compare_same_type(other);
}

ex covdiff::op(size_t i) const{
 if(i==0)
return arg;
if(seq.nops() > 0 && i != 0) 
  return seq[i-1];
 if(seq.nops()==0 && i==1)
return index;
else
throw std::range_error("covdiff::let_op(): no such operand");

}


size_t covdiff::nops() const
{
  if(seq.nops()>0)
    return seq.nops()+1;
  else return 2;
}
ex & covdiff:: let_op(size_t i)
{
  ensure_if_modifiable();
  if(i==0)
return arg;
if(seq.nops() > 0 && i != 0) 
  return seq[i-1];
 if(seq.nops()==0 && i==1)
return index;
else
throw std::range_error("covdiff::let_op(): no such operand");
}





ex covdiff::subs(const exmap & m, unsigned options) const
{
        size_t num = nops();
        if (num) {

                // Substitute in subexpressions
                for (size_t i=0; i<num; i++) {
                        const ex & orig_op = op(i);
                        const ex & subsed_op = orig_op.subs(m, options);
                        if (!are_ex_trivially_equal(orig_op, subsed_op)) {

                                // Something changed, clone the object
                                basic *copy = duplicate();
                                copy->setflag(status_flags::dynallocated);
                                copy->clearflag(status_flags::hash_calculated | 
status_flags::expanded);

                                // Substitute the changed operand
                                copy->let_op(i++) = subsed_op;

                                // Substitute the other operands
                                for (; i<num; i++)
                                        copy->let_op(i) = op(i).subs(m, options);

                                // Perform substitutions on the new object as a
whole
                                return copy->subs_one_level(m, options);
                        }
                }
        }

        // Nothing changed or no subexpressions
        return subs_one_level(m, options);
}

ex Derivative_helper(ex e, ex var)
{
 lst iv;
 if(is_a<covdiff>(e)){
   for(size_t i=1; i<e.nops();i++)
     iv.append(e.op(i));
	       iv.append(var);
	       return covdiff(e.op(0),iv);
 }else
   if(is_a<indexed>(e)){
 if(is_a<su3f>(e.op(0)) && e.nops()==4)
 return 0;
else
   return covdiff(e,var);
   }
   else 
     return 0;
}




struct make_diff : public map_function{
  ex  var;
  make_diff(ex var_): var(var_){}
  ex operator()(const ex & e)
{
if(is_a<add>(e))
      return e.map(*this);
 else
 if(is_a<mul>(e))
	{
	  ex result =0;
	  for(size_t i=0;i<e.nops();++i)
	    {
	      result +=e.subs(e.op(i)==Derivative_helper(e.op(i),var));
	    };
	  return result;
	}else
	  return Derivative_helper(e,var);
}
};
ex D(ex e, ex index){
  ex result;
  make_diff D(index);
   return D(e);
}

-------------------------------------------------
This mail sent through IMP: http://horde.org/imp/