InterpolatedVectorDerivative.cpp

Demonstration of the interpolated derivative of a vector field.

Demonstration of the interpolated derivative of a vector field.

#include <vector/CubicIpol.hpp>
#include <vector/LinearIpol.hpp>
#include <vector/MultiArray.hpp>
#include <vector/Interpolate.hpp>
#include <eagle/Vector.hpp>
#include <eagle/Matrix.hpp>
#include <eagle/QuadraticMatrix.hpp>
 
#include <time.h>
 
using namespace std;
 
using namespace Eagle;
using namespace Fiber;
 
template class LowerTriangular<3, double>;
 
 
// Demonstration on how to compute the partial derivative of
// the nth component of a multi-valued (e.g. vector or tensor)
// field by the mth coordinate function.
static int DemonstrateTensorDerivative()
{
// symbolic names for coordinates
        enum { x=0, y=1, z=2 };
 
// untyped tensorial object 
typedef FixedArray<double,6> uTensor;
 
// typed tensorial object:
//    Symmetric tensor
typedef LowerTriangular<3, double> Tensor;
 
//    General tensor of rank two.
typedef Quadratic<3, double> gTensor;
 
// The data storage of the tensor field is an automatic variable.
Tensor TensorStorage[1000];
 
// This tensor field, consisting of 1000 points, is defined over 
// a three-dimensional base manifold with 10 vertices in each dimension.
MultiArray<3, Tensor> G(TensorStorage, MIndex(10,10,10) );
 
// Initialize all elements through the MultiArray interface
Tensor InitValue;
        InitValue(0,0) = InitValue(1,1) = InitValue(2,2) = 1;
        InitValue(0,1) = InitValue(1,2) = InitValue(2,0) = 0;           
        G.set( InitValue );
 
// Define some point of interesting.
MultiIndex<3> I;
        I = 5,6,7;
 
        // modify the data values there
        G[ I ] (2,2) = 2;
 
        // Determine the vectorial index of the tensorial component
int     Component = Tensor::mem_index(2,2);
 
// Define field referring to the third component of the tensor field.
// Note that we need to view the tensorial data as vector here, thus
// requiring the vec() call. 
Iterator<double> GnmPtr( G.count()*Tensor::SIZE, G.ptr(), Component);
 
        cout << G[I];
 
// Define a multidimensional array on the view pointer.
// Note that the view pointer could also be constructed as local argument.
typedef MultiArray<3, double> Gnm_t;
        Gnm_t Gnm( GnmPtr, G.Size() );
 
        {
        FixedArray<double,3> point;
 
//      Interpolate<3, double, LinearIpol<double> , double, NoDelimiter<double>, x> 
        Interpolate<3, double, CubicIpol<double>  , double, NoDelimiter<double>, x> 
                Gnm_k(Gnm, point);
 
                for(double fr=0; fr<2; fr+=0.1)
                {
                point = 5.1, 6.2, 7 + fr;
                cout << "Cubic interpolated derived tensor component Gmn,k " << point << " " << Gnm_k.eval() <<  endl
                     ;
                }
        }
 
 
        // Starting point
VVector<3, double> q;
        q = 5.1, 6.8, 5.03;
 
        // Velocity
Row<3, double> v; 
        v = 0.01,0.02,.2;
 
clock_t start, finish;
        start = clock();
        for(int is = 0; is<40000; is++)
        {
        //
        // Computation of coordinate acceleration from metric tensor at specific point.
        //
 
        // path parameter increment
double  ds = 0.001;
 
#ifndef EXPLICIT_CHRISTOFFEL_OPERATIONS
//FixedArray<3, FixedArray<3, FixedArray<3, double> > > g_ij_k; 
 
LowerTriangular<3, FixedArray<double,3> > g_ij_k;
 
#if 0
TensorDerivative<MultiArray<3, Tensor, Tensor*> > dG(G, q);
 
        dG.eval( g_ij_k );
#elif 0
 
//const char C[4] = "xyz";
        {for(int j=0; j<3; j++)
         for(int i=j; i<3; i++)
         {
         int    Component = Tensor::mem_index(i,j);
         Gnm_t Gnm( &G.storage().vec(), Component, G.Size() );
#if 1
         Derivator<MultiArray<3, double, ViewPtr<double> >, CubicIpol<double>, double, NoDelimiter<double> > 
                        D(Gnm, q);
 
                D.eval( g_ij_k(i,j) );
#else
#if 1
         Interpolate<3, double, Gnm_t::Storage_t, CubicIpol<double>       , double, NoDelimiter<double>, 0> Gij_0(Gnm, q );
         Interpolate<3, double, Gnm_t::Storage_t, CubicIpol<double>       , double, NoDelimiter<double>, 1> Gij_1(Gnm, q );
         Interpolate<3, double, Gnm_t::Storage_t, CubicIpol<double>       , double, NoDelimiter<double>, 2> Gij_2(Gnm, q );
#else
         Interpolate<3, double, Gnm_t::Storage_t, LinearIpol<double>       , double, NoDelimiter<double>, 0> Gij_0(Gnm, q );
         Interpolate<3, double, Gnm_t::Storage_t, LinearIpol<double>       , double, NoDelimiter<double>, 1> Gij_1(Gnm, q );
         Interpolate<3, double, Gnm_t::Storage_t, LinearIpol<double>       , double, NoDelimiter<double>, 2> Gij_2(Gnm, q );
#endif
#if 0
                g_ij_k[i][j][0] = Gij_0.eval();
                g_ij_k[i][j][1] = Gij_1.eval();
                g_ij_k[i][j][2] = Gij_2.eval();
#else
                g_ij_k(i,j)[0] = Gij_0.eval();
                g_ij_k(i,j)[1] = Gij_1.eval();
                g_ij_k(i,j)[2] = Gij_2.eval();
#endif
 
        //      if (!is)
        //      {
        //              cout << "g_"<< C[i] << C[j] << ",x=" << g_ij_k[i][j][0] << " ";
                //      cout << "g_"<< C[i] << C[j] << ",y=" << g_ij_k[i][j][1] << " ";
                //      cout << "g_"<< C[i] << C[j] << ",z=" << g_ij_k[i][j][2] << " ";
        //      }
#endif
         }
        }
#endif
 
 
//FixedArray<3, FixedArray<3, FixedArray<3, double> > > G_ijk; 
 
LowerTriangular<3, FixedArray<double,3> > G_ijk;
        // Gij_k = gik,j + gjk,i - gij,k
 
        {for(int j=0; j<3; j++)
         for(int i=j; i<3; i++)
         for(int k=0; k<3; k++)
         {
//               G_ijk[i][j][k] = g_ij_k[i][k][j] + g_ij_k[j][k][i] - g_ij_k[i][j][k];
                 G_ijk(i,j)[k] = g_ij_k(i,k)[j] + g_ij_k(j,k)[i] - g_ij_k(i,j)[k];
         }
        }
 
Interpolate<3, Tensor, CubicIpol<Tensor>, double> Gfield(G, q );
 
 
const Tensor&Glocal = Gfield.eval();
 
Column<3, double>       dv;
 
        // Compute coordinate co-acceleration via Lower Christoffels
        //
        // dv_k = G_{ij_k} v^i v^j
        //
        {for(int k=0; k<3; k++)
        {
                dv[k] = 0;
                for(int j=0; j<3; j++)
                {
                double  G_ijk_vi = 0;
                        for(int i=0; i<3; i++)
                        {
                                G_ijk_vi += G_ijk(i,j)[k] * v[i];
                        }
                        G_ijk_vi *= v[j];
                        dv[k] += G_ijk_vi; 
                }
        }
        }
 
Quadratic<3,double> Glocal_q = Glocal;
Quadratic<3,double>::PermutationVector_t PV;
        Glocal_q.GaussDecompose(PV); 
 
        // coordinate acceleration, \ddot q
Column<3, double> a;
 
        Glocal_q.GaussSolve(PV, dv, a );
 
        cout << " Coordinate Acceleration:" << a << endl;
 
#else // EXPLICIT_CHRISTOFFEL_OPERATIONS
 
ChristoffelField<MultiArray<3, Tensor, Tensor*> > Chris(G,q);
        Chris.getAcceleration(a, v);
#endif
 
        // v == dq/ds
        // a == dv/ds; 
 
        #pragma message "this next operator should work - CHECK!"
        //q += v*ds; 
 
        v += a*ds;
 
                if (is %1000 == 0)
                        cout << "v="<<v<<"\t q="<<q<<endl;
 
        }
        finish = clock();
        printf( "Geodesic computation: %2.1f seconds\n", (double)(finish - start) / CLOCKS_PER_SEC );
        
        cout << "Final: v="<<v<<"\t q="<<q<<endl;
        cout << "EXPECT v=(0.00628477,0.088151,0.180334)  q=(5.3987,9.70163,12.4076)" << endl;
 
/*
----- P4 1.6GHz, Cubic Christoffel Ipol, Intel C++, Optimize:
 
Geodesic computation: 5.7 seconds
Final: v=(0.00628477,0.088151,0.180334)  q=(5.3987,9.70163,12.4076)
 
----> Employing symmetry on Christoffel:
Geodesic computation: 4.0 seconds
Final: v=(0.00628477,0.088151,0.180334)  q=(5.3987,9.70163,12.4076)
 
----- P4 1.6GHz, Linear Christoffel Ipol, Intel C++, Optimize:
Geodesic computation: 1.8 seconds
Final: v=(0.00160945,0.00741571,0.0737203)       q=(5.38425,7.92392,10.8191)
 
----- P4 1.6GHz, Cubic Christoffel Ipol, Intel C++, Debug:
 
Geodesic computation: 56.5 seconds
Final: v=(0.00628477,0.088151,0.180334)  q=(5.3987,9.70163,12.4076)
        */ 
        return 1;
}
 
 
// Demonstration the computation of geodesics in a numerical tensor field
// and the numerical precision via back-integration to the starting point.
static int DemonstrateGeodesics()
{
// symbolic names for coordinates
        enum { x=0, y=1, z=2 };
 
// typed tensorial object:
typedef LowerTriangular<3, double> Tensor;
 
// The data storage of the tensor field is an automatic variable.
Tensor TensorStorage[1000];
 
// This tensor field, consisting of 1000 points, is defined over 
// a three-dimensional base manifold with 10 vertices in each dimension.
MultiArray<3, Tensor> G(TensorStorage, MIndex(10,10,10) );
 
// Initialize all elements through the MultiArray interface
Tensor  InitValue;
        InitValue(0,0) = InitValue(1,1) = InitValue(2,2) = 1;
        InitValue(0,1) = InitValue(1,2) = InitValue(2,0) = 0;           
        G.set( InitValue );
 
        // Modify the data values at some point of interest
        G[ MIndex(5,6,7) ] (2,2) = 2;
 
 
 
        // Starting point
VVector<3, double> q0, q;
        q = 5.1, 6.8, 5.03;
        q0 = q;
 
        // Velocity
Row<3, double> v0, v; 
        v = 0.01,0.02,.2;
        v0 = v;
 
#if 0
ChristoffelField<MultiArray<3, Tensor, Tensor*> > Chris(G,q); 
 
// EXCERPT: 
#if 1
 {
 typedef Tensor metric_type; 
 typedef Tensor* MetricStorage_t; 
 typedef CubicIpol<metric_type> MetricInterpol; 
 typedef double CoordinateType;
  FixedArray<3, CoordinateType> p;
  Interpolate<3, metric_type, MetricStorage_t, MetricInterpol, CoordinateType> MetricField(G, p); 
 
  NoDelimiter<metric_type> ND; 
 
//  metric_type m = G[0]; 
  metric_type m = MetricField[0];
 
//  metric_type m = ND.limit( MetricField[0] ); 
 
//  metric_type m = MetricField.myDelimiter.limit( MetricField[0] ); 
 
//MetricInterpol::interpolate( MetricField, p[2],  G.Size[2], MetricField.myDelimiter ); 
 
//IpolDerivative<metric_type, MetricInterpol, false>::interpolate(MetricField, p[2],  G.Size[2], MetricField.myDelimiter ); 
 
//      metric_type g = MetricField.eval();
 }
 
#else
 
clock_t start, finish;
        start = clock();
        {for(int is = 0; is<40000; is++)
        {
                // path parameter increment
        double  ds = 0.001;
 
                // coordinate acceleration, \ddot q
        Row<3, double> a;
 
                Chris.getCoordinateAcceleration(a, v);
 
                // Euler step
                q += v*ds;
                v += a*ds;
 
                if (is %1000 == 0)
                        cout << "v="<<v<<"\t q="<<q<<endl;
        }}
        finish = clock();
        printf( "Geodesic computation: %2.1f seconds\n", (double)(finish - start) / CLOCKS_PER_SEC );
 
        {for(int is = 0; is<40000; is++)
        {
                // path parameter increment
        double  ds = -0.001;
 
                // coordinate acceleration, \ddot q
        Row<3, double> a;
 
                Chris.getCoordinateAcceleration(a, v);
 
                // Euler step
                q += v*ds;
                v += a*ds;
 
                if (is %1000 == 0)
                        cout << "v="<<v<<"\t q="<<q<<endl;
        }}
 
        cout << "q=" << q << "\t q0=" << q0 << " Diff: " << q-q0 << endl;
        cout << "v=" << v << "\t v0=" << v0 << " Diff: " << v-v0 <<  endl;
 
double  Qerror = norm2(q-q0);
        cout << "Qerror =" << Qerror  << endl;
        assert( Qerror < 1E-8 );
#endif
#endif 
        return 1;
}
 
 
int dVInterpolCheck()
{
int r=0; 
//      r+=DemonstrateTensorDerivative(); 
//      r+=DemonstrateGeodesics(); 
        return r;
}