Geodesics .hpp

#ifndef __GEODESCICS_HPP_
#define __GEODESCICS_HPP_

#include <fish/fiber/vector/CubicIpol.hpp>
#include <eagle/QuadraticMatrix.hpp> 
#include <eagle/Coordinates.hpp>
#include <fish/fiber/baseop/LocalFromWorldPoint.hpp>
#include <ode/dop853.hpp>

namespace FieldLines
{

struct Geodesic{};


template<class T>
class ChristoffelXYZ : public MemCore::ReferenceBase< ChristoffelXYZ<T> >
{
        typedef typename T::Chart_t Chart_t;
        typedef typename Eagle::Coordinates<Chart_t>::point point;
        // store metric tensor
        T g;
        
        Eagle::Quadratic<T::Dims, double> g_inv;
 
        // store precomputed finite differences of g // array size 3
        T g_diff[T::Dims];


        bool ok;
        // consts for access
//      enum { MU = 0, NU, LAMBDA  };
//      enum { X = 0, Y, Z};

public:
        ChristoffelXYZ( const T&gP, const T&g_xP, const T&g_yP, const T&g_zP,
                        const double dxP, const double dyP, const double dzP  )
        : MemCore::ReferenceBase< ChristoffelXYZ >(this)
        , g(gP)
        , ok(true)
        {
                if(T::Dims == 3)
                {
                        enum { X = 0, Y, Z};

                if(dxP != 0.0)
                        for(unsigned i = 0; i < T::SIZE; i++)
                                g_diff[X][i] = ( g_xP[i] - g[i] ) / dxP; 

                if(dyP != 0.0)
                        for(unsigned i = 0; i < T::SIZE; i++)
                                g_diff[Y][i] = ( g_yP[i] - g[i] ) / dyP; 

                if(dzP != 0.0)
                        for(unsigned i = 0; i < T::SIZE; i++)
                                g_diff[Z][i] = ( g_zP[i] - g[i] ) / dzP;

                } 

                if(T::Dims == 4)
                {
                        enum {  t = Eagle::STA::CartesianChart4D<>::T, 
                                X = Eagle::STA::CartesianChart4D<>::X,
                                Y = Eagle::STA::CartesianChart4D<>::Y,
                                Z = Eagle::STA::CartesianChart4D<>::Z }; 

                        if(dxP != 0.0)
                                for(unsigned i = 0; i < T::SIZE; i++)
                                        g_diff[X][i] = ( g_xP[i] - g[i] ) / dxP; 

                        if(dyP != 0.0)
                                for(unsigned i = 0; i < T::SIZE; i++)
                                        g_diff[Y][i] = ( g_yP[i] - g[i] ) / dyP; 

                        if(dzP != 0.0)
                                for(unsigned i = 0; i < T::SIZE; i++)
                                        g_diff[Z][i] = ( g_zP[i] - g[i] ) / dzP;
                        
                }
        }

        ChristoffelXYZ( const T&gP, const T&g_dxP, const T&g_dyP, const T&g_dzP )
        : MemCore::ReferenceBase< ChristoffelXYZ<T> >(this)
        , g(gP)
        , ok(true)      
        {
//              std::cout << "hello?";
                if(T::Dims == 3)
                {
                        enum { X = 0, Y, Z};
                        g_diff[X] = g_dxP; 
                        g_diff[Y] = g_dyP; 
                        g_diff[Z] = g_dzP; 

                        try
                        {
                                Eagle::Quadratic<T::Dims, double> q (g); 

                                Eagle::ComputeInverse(g_inv, q); 

                                // Eagle::Matrix<T::Dims,T::Dims, double> &mg = q ; 
                                // Eagle::Matrix<T::Dims,T::Dims, double> &mg_i = g_inv; 
                                // Eagle::Matrix<T::Dims,T::Dims, double> I = mg * mg_i; 
                                // std::cout << "controll inversion: " << I << std::endl; 
                        } 
                        catch(const Eagle::DegeneratedMatrix&)
                        {
                                puts("ChristoffelXYZ<>::ChristoffelXYZ ERROR Degenerated Matrix in inversion."); fflush(stdout); 
                                ok = false;
                        } 
                        
                } 
                if(T::Dims == 4)
                {
                        enum {  t = Eagle::STA::CartesianChart4D<>::T, 
                                X = Eagle::STA::CartesianChart4D<>::X,
                                Y = Eagle::STA::CartesianChart4D<>::Y,
                                Z = Eagle::STA::CartesianChart4D<>::Z }; 

                        g_diff[t].set(0.0);
                        g_diff[X] = g_dxP; 
                        g_diff[Y] = g_dyP; 
                        g_diff[Z] = g_dzP; 

                        g_inv(g);
                } 



        }

        bool isOk(){ return ok;}

        void  getMetricfromQuadratic(const Eagle::Quadratic<T::Dims,double>&q, T&g)
                {
                        for(unsigned i = 0; T::SIZE; i++)
                                g[i] = q[i];
                }



        // christoffel symbols function for on demand symbol computing, later there may be the need to 
        // store them and reuse them when possible.
        const double operator()(const unsigned mu, const unsigned nu, const unsigned lambda) const 
        {
        double gamma = 0.0; 

                if(mu > T::Dims-1 || nu > T::Dims-1 || lambda > T::Dims-1)
                        return 0.0; 

                for(unsigned i = 0; i < T::Dims; i++)
                        gamma += g_inv(lambda, i) * ( g_diff[nu](mu,i) + g_diff[mu](nu,mu) - g_diff[i](mu, nu) ); 

                gamma /= 0.5;

                return gamma;
        }

        const T& getG() { return g; }
        
        void speak() const
       {
               if(T::Dims == 3)
                       std::cout << "Christoffel content:\n" << "g: " << g << "\ng': " << g_diff[0] << g_diff[1] << g_diff[2]; 
               if(T::Dims == 4)
                       std::cout << "Christoffel content:\n" << "g: " << g << "\ng': " << g_diff[0] << g_diff[1] << g_diff[2]<< g_diff[3];

       }

};

template<class T>
typename Eagle::Coordinates<typename T::Chart_t >::vector getQddot( RefPtr< ChristoffelXYZ<T> >& Gamma, 
                                                                    typename Eagle::Coordinates<typename T::Chart_t >::vector&q_dot )
{
typename Eagle::Coordinates<typename T::Chart_t >::vector q_ddot; 

        q_ddot.set(0.0);

        // q-ddot^lamba = - Gamma^lambda_mu_nu * q-dot^mu * q-dot^nu 

        for(index_t lambda = 0; lambda < T::Dims; lambda++)
                for( index_t mu = 0; mu < T::Dims; mu++)
                        for(index_t nu = 0; nu < T::Dims; nu++)
                                q_ddot[lambda] -= (*Gamma)( mu, nu, lambda ) * q_dot[mu] * q_dot[nu]; 
        return q_ddot;
}

template<class T>
bool getChristoffelXYZ( const typename Eagle::Coordinates<typename T::Chart_t >::point& CurrentPoint, 
                        const MemCore::RefPtr<LocalFromWorldPoint>&LocalPointFinder, 
                        const FieldSelector&FieldSelection, const double step_size, 
                        RefPtr< ChristoffelXYZ<T> >& Christoffels,
                        T& CurrentMetric,
                        pair<Eagle::PhysicalSpace::point, string>& localPoint ) 
        {
        typedef typename T::Chart_t Chart_t;
        typedef typename Coordinates<Chart_t>::point point; 

                Eagle::PhysicalSpace::point current (CurrentPoint[0+T::Dims-3], CurrentPoint[1+T::Dims-3], CurrentPoint[2+T::Dims-3]); 

        bool test = LocalPointFinder->get( current, localPoint ); 

                if(!test)
                {
                        puts("getChristoffelXYZ()1 error finding local point at:"); 
                        cout << current << endl;
                        return false;
                } 
                
        const string&FragName = localPoint.second; 
        const Eagle::PhysicalSpace::point &float_index = localPoint.first; 

//              cout << "FieldLines::getChristoffelXYZ local point: " << float_index; 
//              cout << "in Fragment:" << FragName << endl;
               
        RefPtr<Field> Field; 

                Field = FieldSelection.getField(); 

//              cout << "getting Field: " << FieldSelection.getFieldName() << endl;
                
                if(!Field) 
                { 
                        //puts("AtomicIntegrator<metric33, Geodesic>doit2(); ERROR: No Field0 found in Time iteration"); 
                        return false; 
                } 


         RefPtr<FragmentID> FragID;

                if(FragName.compare("") != 0 )
                        FragID   = new FragmentID( FragName ); 

                RefPtr<CreativeArrayBase> cBase = Field->getCreator( FragID ); 
                assert(cBase); 

                RefPtr<MemBase> mBase = cBase->create(); 
                assert(mBase); 

//              mBase.speak("bla");

                RefPtr<MemArray<3, T> > ToIntArr0 = mBase; 
                ToIntArr0   = Field->getCreator( FragID )->create(); 

                assert( ToIntArr0 );



                        // NOTE: Here optionally use another interpolator, eventually cubic, or 
                        //       self-aligning for eigenvector fields.
                        // 

//              cout << "before interpol" << endl;

        T InterpolData0_x, InterpolData0_y, InterpolData0_z; 

        Interpolate<3, T, CubicIpol<T> > myField0  ( *ToIntArr0,   float_index); 
        Interpolate<3, T, CubicIpol<T>, double, NoDelimiter<T>, 0  > myField0_x( *ToIntArr0, float_index); 
        Interpolate<3, T, CubicIpol<T>, double, NoDelimiter<T>, 1  > myField0_y( *ToIntArr0, float_index); 
        Interpolate<3, T, CubicIpol<T>, double, NoDelimiter<T>, 2  > myField0_z( *ToIntArr0, float_index); 
                  CurrentMetric   = myField0.eval(); 
                  InterpolData0_x = myField0_x.eval(); 
                  InterpolData0_y = myField0_y.eval(); 
                  InterpolData0_z = myField0_z.eval(); 

//                cout << "-----------------------" << endl << CurrentMetric << InterpolData0_x << InterpolData0_y << InterpolData0_z; 

                  Christoffels = new ChristoffelXYZ<T>( CurrentMetric, InterpolData0_x, InterpolData0_y, InterpolData0_z ); 

                  if( Christoffels->isOk() == false)
                  {
                          puts("getChristoffelXYZ()<T>:warning, christoffels probably degenerated matrix."); 
                          return false;
                  }

                  return true;
        }


template<class T>
struct  GeodesicDifferentialEquation : public Traum::VirtualSecondOrderDiffEquation <double>
{
typedef typename T::Chart_t Chart_t;
typedef typename Eagle::Coordinates<Chart_t>::point point;
typedef typename Eagle::Coordinates<Chart_t>::vector tvector;

        
//typedef ::Eagle::metric33 metric33;

RefPtr<LocalFromWorldPoint> LocalPointFinder;   
FieldSelector FieldSelection;

//MultiIndex<3>       FieldSize;
        
//FixedArray<double, 3> DomainStart,
//                    DomainEnd;

double step_size; // define stepsize for finite differential derivation of the metric tensors used for ChristoffelsXYZ33
double scale;

bool valid_point;

        pair<Eagle::PhysicalSpace::point, std::string>localPoint;

T CurrentMetric;
        Eagle::tvector3 DataDir;
RefPtr<ChristoffelXYZ<T> > Gamma;

        enum { N = 6 };
        typedef double real;
        
        GeodesicDifferentialEquation()
        {}
        
        virtual ~GeodesicDifferentialEquation() {}
        
        struct RealArray_t : FixedArray<real, N>               
        {
                void init(int) {}
        };
        
        struct DOPVarsArray_t : FixedArray<Traum::dop_vars<real>, N >
        {
                void init(int) {}
        };
        
        void Accel(real s, const real *x, const real *v, real *d2x_ds2)
        {
                point Q; 
                tvector q_dot; 

                        for(int k=0; k < point::Dims; k++)
                        {
                                Q[k] = x[k]; 
                                q_dot[k] = v[k];
                        } 
//              printf("Q: %f, %f, %f\n", Q[0], Q[1], Q[2] ); 
//              printf("V: %f, %f, %f\n", v[0], v[1], v[2] ); 

        bool test = FieldLines::getChristoffelXYZ<T>( Q, LocalPointFinder, FieldSelection, step_size, 
                                                     Gamma,
                                                     CurrentMetric,
                                                     localPoint ); 

                if(!test)
                {
                        //puts("GeodesicDifferentialEquation::DifferentialEquation::Accel() No valid point found!");fflush(stdout); 
                        return;
                } 

                tvector qdd = getQddot( Gamma, q_dot); 

                for(int k=0; k < point::Dims; k++)
                        d2x_ds2[k] = qdd[k];

        }

};

// double linearInterpolate(const double t0, const double d0, const double t1, const double d1, const double t)
// {
//      return d0 + (d1 - d0) / (t1 - t0) * (t - t0);
// }

// metric33 linearInterpolate(const double t0, const metric33&m0, const double t1, const metric33&m1, const double t )
// {
// metric33 tmp; 
//      for(unsigned i = 0; i < metric33::SIZE; i++)
//              tmp[i] = linearInterpolate(t0, m0[i], t1, m1[i], t); 

//      return tmp;
// }


}
#endif // __GEODESCICS_HPP_