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 <vector/Vector.hpp> #include <vector/Matrix.hpp> #include <vector/QuadraticMatrix.hpp> #include <time.h> using namespace std; 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<6, double> 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<3, double> 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<3, double> > 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<3, double> > 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; 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; }
1.6.2