00001 #ifndef __GEODESCICS_HPP_
00002 #define __GEODESCICS_HPP_
00003
00004 #include <fish/fiber/vector/CubicIpol.hpp>
00005 #include <fish/fiber/vector/LinearIpol.hpp>
00006 #include <eagle/QuadraticMatrix.hpp>
00007 #include <eagle/Coordinates.hpp>
00008 #include <fish/fiber/baseop/LocalFromWorldPoint.hpp>
00009 #include <ode/dop853.hpp>
00010
00011
00012
00013 namespace Fiber
00014 {
00015
00016 template <>
00017 struct LinearIpolZeroDerivativeTrait<metric33>
00018 {
00019 static metric33 zero()
00020 {
00021 metric33 m;
00022 m.set(0.0);
00023 return m;
00024 }
00025 };
00026
00027 template <>
00028 struct LinearIpolZeroDerivativeTrait<metric44>
00029 {
00030 static metric44 zero()
00031 {
00032 metric44 m;
00033 m.set(0.0);
00034 return m;
00035 }
00036 };
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060 }
00061
00062
00063 namespace FieldLines
00064 {
00065
00066 struct Geodesic{};
00067
00068 struct GeodesicSpatial{};
00069
00070 template<class T>
00071 double getDet(const Eagle::Quadratic<T::Dims,double>&q)
00072 {
00073 return Eagle::Determinantor<T::Dims>::compute(q);
00074 }
00075
00076 template<class T>
00077 class ChristoffelXYZ : public MemCore::ReferenceBase< ChristoffelXYZ<T> >
00078 {
00079 typedef typename T::Chart_t Chart_t;
00080 typedef typename Eagle::Coordinates<Chart_t>::point point;
00081
00082 T g;
00083
00084 Eagle::Quadratic<T::Dims, double> g_inv;
00085
00086
00087 T g_diff[T::Dims];
00088
00089
00090 bool ok;
00091
00092
00093
00094
00095 public:
00096 ChristoffelXYZ( const T&gP, const T&g_xP, const T&g_yP, const T&g_zP,
00097 const double dxP, const double dyP, const double dzP )
00098 : MemCore::ReferenceBase< ChristoffelXYZ >(this)
00099 , g(gP)
00100 , ok(true)
00101 {
00102 std::cout << "aber jetzt nit im ersnst!!!";
00103 if(T::Dims == 3)
00104 {
00105 enum { X = 0, Y, Z};
00106
00107 if(dxP != 0.0)
00108 for(unsigned i = 0; i < T::SIZE; i++)
00109 g_diff[X][i] = ( g_xP[i] - g[i] ) / dxP;
00110
00111 if(dyP != 0.0)
00112 for(unsigned i = 0; i < T::SIZE; i++)
00113 g_diff[Y][i] = ( g_yP[i] - g[i] ) / dyP;
00114
00115 if(dzP != 0.0)
00116 for(unsigned i = 0; i < T::SIZE; i++)
00117 g_diff[Z][i] = ( g_zP[i] - g[i] ) / dzP;
00118
00119 }
00120
00121 if(T::Dims == 4)
00122 {
00123 enum { t = Eagle::STA::CartesianChart4D<>::T,
00124 X = Eagle::STA::CartesianChart4D<>::X,
00125 Y = Eagle::STA::CartesianChart4D<>::Y,
00126 Z = Eagle::STA::CartesianChart4D<>::Z };
00127
00128 if(dxP != 0.0)
00129 for(unsigned i = 0; i < T::SIZE; i++)
00130 g_diff[X][i] = ( g_xP[i] - g[i] ) / dxP;
00131
00132 if(dyP != 0.0)
00133 for(unsigned i = 0; i < T::SIZE; i++)
00134 g_diff[Y][i] = ( g_yP[i] - g[i] ) / dyP;
00135
00136 if(dzP != 0.0)
00137 for(unsigned i = 0; i < T::SIZE; i++)
00138 g_diff[Z][i] = ( g_zP[i] - g[i] ) / dzP;
00139
00140 }
00141 }
00142
00143 ChristoffelXYZ( const T&gP, const T&g_dxP, const T&g_dyP, const T&g_dzP )
00144 : MemCore::ReferenceBase< ChristoffelXYZ<T> >(this)
00145 , g(gP)
00146 , ok(true)
00147 {
00148
00149
00150 if(T::Dims == 3)
00151 {
00152 enum { X = 0, Y, Z};
00153 g_diff[X] = g_dxP;
00154 g_diff[Y] = g_dyP;
00155 g_diff[Z] = g_dzP;
00156 }
00157 if(T::Dims == 4)
00158 {
00159 enum { t = Eagle::STA::CartesianChart4D<>::T,
00160 X = Eagle::STA::CartesianChart4D<>::X,
00161 Y = Eagle::STA::CartesianChart4D<>::Y,
00162 Z = Eagle::STA::CartesianChart4D<>::Z };
00163
00164 g_diff[t].set(0.0);
00165 g_diff[X] = g_dxP;
00166 g_diff[Y] = g_dyP;
00167 g_diff[Z] = g_dzP;
00168 }
00169 try
00170 {
00171 Eagle::Quadratic<T::Dims, double> q (g);
00172
00173 Eagle::ComputeInverse(g_inv, q);
00174
00175 Eagle::Matrix<T::Dims,T::Dims, double> &mg = q ;
00176 Eagle::Matrix<T::Dims,T::Dims, double> &mg_i = g_inv;
00177 Eagle::Matrix<T::Dims,T::Dims, double> I = mg * mg_i;
00178
00179 }
00180 catch(const Eagle::DegeneratedMatrix&)
00181 {
00182 puts("ChristoffelXYZ<>::ChristoffelXYZ ERROR Degenerated Matrix in inversion."); fflush(stdout);
00183 ok = false;
00184 }
00185 }
00186
00187 bool isOk(){ return ok;}
00188
00189 void getMetricfromQuadratic(const Eagle::Quadratic<T::Dims,double>&q, T&g)
00190 {
00191 for(unsigned i = 0; T::SIZE; i++)
00192 g[i] = q[i];
00193 }
00194
00195
00196
00197
00198
00199 const double operator()(const unsigned mu, const unsigned nu, const unsigned lambda) const
00200 {
00201 double gamma = 0.0;
00202
00203 if(mu > T::Dims-1 || nu > T::Dims-1 || lambda > T::Dims-1)
00204 return 0.0;
00205
00206 for(unsigned i = 0; i < T::Dims; i++)
00207 gamma += g_inv(lambda, i) * ( g_diff[nu](mu,i) + g_diff[mu](nu,i) - g_diff[i](mu, nu) );
00208
00209 gamma /= 2;
00210
00211 return gamma;
00212 }
00213
00214 const double gammaDebug(const unsigned mu, const unsigned nu, const unsigned lambda) const
00215 {
00216 double gamma = 0.0;
00217
00218 if(mu > T::Dims-1 || nu > T::Dims-1 || lambda > T::Dims-1)
00219 return 0.0;
00220 printf( "mu nu lam %d %d %d", mu, nu, lambda);
00221 for(unsigned i = 0; i < T::Dims; i++)
00222 {
00223
00224 double tmp = g_inv(lambda, i) * ( g_diff[nu](mu,i) + g_diff[mu](nu,i) - g_diff[i](mu, nu) );
00225 gamma += tmp;
00226 printf( "g_inv: %f, g1: %f, g2:% f, g3: %f, term: %f\n", g_inv(lambda, i), g_diff[nu](mu,i), g_diff[mu](nu,i), g_diff[i](mu, nu), tmp);
00227 }
00228 gamma /= 2;
00229 printf( "gamma: %f\n", gamma);
00230 return gamma;
00231 }
00232
00233
00234 const T& getG() { return g; }
00235
00236 void speak() const
00237 {
00238
00239 std::cout << "Christoffel content:\n" << "g " << endl << g << "g_inv: " << endl << g_inv;
00240
00241 if(T::Dims == 3)
00242 std::cout << "g,x" << std::endl <<g_diff[0]
00243 << "g,y" << std::endl <<g_diff[1]
00244 << "g,z" << std::endl <<g_diff[2];
00245
00246 if(T::Dims == 4)
00247 std::cout << "g,t" << std::endl <<g_diff[0]
00248 << "g,x" << std::endl <<g_diff[1]
00249 << "g,y" << std::endl <<g_diff[2]
00250 << "g,z" << std::endl <<g_diff[3];
00251
00252 cout << "Gamma (T)|X|Y|Z|:" << endl;
00253
00254 for( unsigned lam = 0; lam < T::Dims; lam++)
00255 {
00256 for( unsigned mu = 0; mu < T::Dims; mu++)
00257 {
00258 for( unsigned nu = 0; nu < T::Dims; nu++)
00259 {
00260 printf("%f ", (*this)(mu,nu,lam));
00261 }
00262 printf("\n");
00263 }
00264 printf("-----------------------------\n");
00265 }
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280 }
00281
00282 };
00283
00284
00287 template<class T>
00288 typename Eagle::Coordinates<typename T::Chart_t >::vector getQddot( RefPtr< ChristoffelXYZ<T> >& Gamma,
00289 typename Eagle::Coordinates<typename T::Chart_t >::vector&q_dot )
00290 {
00291 typename Eagle::Coordinates<typename T::Chart_t >::vector q_ddot;
00292
00293 q_ddot.set(0.0);
00294
00295
00296
00297
00298
00299 for(index_t lambda = 0; lambda < T::Dims; lambda++)
00300 for( index_t mu = 0; mu < T::Dims; mu++)
00301 for(index_t nu = 0; nu < T::Dims; nu++)
00302 q_ddot[lambda] -= (*Gamma)( mu, nu, lambda ) * q_dot[mu] * q_dot[nu];
00303
00304
00305
00306 return q_ddot;
00307 }
00308
00309
00310 template<class TT>
00311 typename Eagle::Coordinates<typename TT::Chart_t >::vector getQdot( const typename Eagle::Coordinates<typename TT::Chart_t >::vector& q_dot,
00312 const TT&g, const bool big = true)
00313 {
00314 if (TT::Dims==3)
00315 return q_dot;
00316
00317 enum { T = STA::CartesianChart4D<>::T,
00318 X = STA::CartesianChart4D<>::X,
00319 Y = STA::CartesianChart4D<>::Y,
00320 Z = STA::CartesianChart4D<>::Z };
00321
00322 const double eps = 1.0e-4;
00323
00324
00325
00326 typename Eagle::Coordinates<typename TT::Chart_t >::vector v = q_dot;
00327
00328 if( !isnormal(v[X]) ) v[X] = 0.0;
00329 if( !isnormal(v[Y]) ) v[Y] = 0.0;
00330 if( !isnormal(v[Z]) ) v[Z] = 0.0;
00331
00332
00333 if(v[X] < eps && v[X] > -eps &&
00334 v[Y] < eps && v[Y] > -eps &&
00335 v[Z] < eps && v[Z] > -eps)
00336 {
00337
00338 v[T] = 1.0;
00339 return v;
00340 }
00341
00342
00343
00344 double A = 2 * ( g(X,T)*v[X] + g(Y,T)*v[Y] + g(Z,T)*v[Z] );
00345 double B = g(X,X)*v[X]*v[X] + g(Y,Y)*v[Y]*v[Y] + g(Z,Z)*v[Z]*v[Z]+
00346 2*( g(X,Y)*v[X]*v[Y] + g(X,Z)*v[X]*v[Z] + g(Y,Z)*v[Y]*v[Z] );
00347
00348 if(g(T,T) == 0.0)
00349 {
00350 puts("getQdot::Error g(T,T) == 0! Division by 0.");
00351 return v;
00352 }
00353
00354 double p = A/g(T,T);
00355 double q= B/g(T,T);
00356
00357 double w = p*p/4 - q;
00358
00359 if(w < 0.0)
00360 {
00361 puts("getQdot::Error direction initialization srqt < 0!!!!");
00362 return v;
00363 }
00364
00365 double sw = sqrt( w );
00366
00367 double v1 = -p/2 + sw;
00368 double v2 = -p/2 - sw;
00369
00370 if(big)
00371 v[T] = (v1 > v2)? v1 : v2;
00372 else
00373 v[T] = (v1 > v2)? v2 : v1;
00374
00375 return v;
00376 }
00377
00378
00381 template<class T>
00382 bool getChristoffelXYZ( const typename Eagle::Coordinates<typename T::Chart_t >::point& CurrentPoint,
00383 const MemCore::RefPtr<LocalFromWorldPoint>&LocalPointFinder,
00384 const FieldSelector&FieldSelection, const double step_size,
00385 RefPtr< ChristoffelXYZ<T> >& Christoffels,
00386 T& CurrentMetric,
00387 pair<Eagle::PhysicalSpace::point, string>& localPoint )
00388 {
00389 typedef typename T::Chart_t Chart_t;
00390 typedef typename Coordinates<Chart_t>::point point;
00391
00392 Eagle::PhysicalSpace::point current (CurrentPoint[0+T::Dims-3], CurrentPoint[1+T::Dims-3], CurrentPoint[2+T::Dims-3]);
00393
00394 bool test = LocalPointFinder->get( current, localPoint );
00395
00396 if(!test)
00397 {
00398 puts("getChristoffelXYZ()1 error finding local point at:");
00399 cout << current << endl;
00400 return false;
00401 }
00402
00403 const string&FragName = localPoint.second;
00404 const Eagle::PhysicalSpace::point &float_index = localPoint.first;
00405
00406
00407
00408 RefPtr<Field> Field;
00409 Field = FieldSelection.getField();
00410
00411 #ifdef VERBOSE
00412 cout << "FieldLines::getChristoffelXYZ local point: " << float_index;
00413 cout << "in Fragment:" << FragName << endl;
00414 cout << "getting Field: " << FieldSelection.getFieldName() << endl;
00415 #endif
00416 if(!Field)
00417 {
00418 puts("AtomicIntegrator<metric33, Geodesic>doit2(); ERROR: No Field0 found in Time iteration");
00419 return false;
00420 }
00421
00422
00423 RefPtr<FragmentID> FragID;
00424
00425 if(FragName.compare("") != 0 )
00426 FragID = new FragmentID( FragName );
00427
00428 RefPtr<CreativeArrayBase> cBase = Field->getCreator( FragID );
00429 assert(cBase);
00430
00431 RefPtr<MemBase> mBase = cBase->create();
00432 assert(mBase);
00433
00434
00435
00436 RefPtr<MemArray<3, T> > ToIntArr0 = mBase;
00437 ToIntArr0 = Field->getCreator( FragID )->create();
00438
00439 if(!ToIntArr0)
00440 {
00441 cout << "getChristoffelXYZ<T> dims:" << T::Dims << endl;
00442 mBase.speak("getChristoffelXYZ<T> ------>");
00443
00444 puts("getChristoffelXYZ()<T>:ERROR retrieving matching field. Maybe Dims of instantiated metric and Dims of field data does not match. 33 44?");
00445 return false;
00446 }
00447
00448 assert( ToIntArr0 );
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458 T g_x, g_y, g_z;
00459 T g_xl, g_xp, g_xxp, g_yp, g_zp, g_yzp, g_yyzp;
00460
00461 Interpolate<3, T, CubicIpol<T> > myField0 ( *ToIntArr0, float_index);
00462 Interpolate<3, T, CubicIpol<T>, double, NoDelimiter<T>, 0 > myField0_x( *ToIntArr0, float_index);
00463 Interpolate<3, T, CubicIpol<T>, double, NoDelimiter<T>, 1 > myField0_y( *ToIntArr0, float_index);
00464 Interpolate<3, T, CubicIpol<T>, double, NoDelimiter<T>, 2 > myField0_z( *ToIntArr0, float_index);
00465
00466
00467
00468
00469
00470 CurrentMetric = myField0.eval();
00471 g_x = myField0_x.eval();
00472 g_y = myField0_y.eval();
00473 g_z = myField0_z.eval();
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502 RefPtr<Fiber::Field> field = FieldSelection.getCartesianPositions();
00503 RefPtr<Fiber::DirectProductMemArray<Eagle::point3,3> > FData = field->getData();
00504 Eagle::point3 delta = FData->Delta();
00505
00506
00507 Eagle::Quadratic<T::Dims, double> q (CurrentMetric);
00508 double det = getDet<T>(q);
00509
00510
00511
00512 if(det < 0.0 && T::Dims == 3)
00513 {
00514 cout << "g:" << q << endl;
00515 cout << "det: " << det << endl;
00516 puts("getChristoffelXYZ()<T>:warning, Det of G < 0.");
00517 return false;
00518 }
00519
00520 g_x /= delta[0];
00521 g_xl /= delta[0];
00522 g_y /= delta[1];
00523 g_z /= delta[2];
00524
00525 #ifdef VERBOSE
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543 #endif
00544 Christoffels = new ChristoffelXYZ<T>( CurrentMetric, g_x, g_y, g_z );
00545
00546 #ifdef VERBOSE
00547 std::cout << "Hilfe Christoffel wo bist du?" << std::endl;
00548 cout << "current point: " << CurrentPoint;
00549 Christoffels->speak();
00550 #endif
00551 if( Christoffels->isOk() == false)
00552 {
00553 puts("getChristoffelXYZ()<T>:warning, christoffels probably degenerated matrix.");
00554 return false;
00555 }
00556
00557 return true;
00558 }
00559
00560 template<class T>
00561 struct GeodesicDifferentialEquation : public Traum::VirtualSecondOrderDiffEquation <double>
00562 {
00563 typedef typename T::Chart_t Chart_t;
00564 typedef typename Eagle::Coordinates<Chart_t>::point point;
00565 typedef typename Eagle::Coordinates<Chart_t>::vector tvector;
00566
00567 RefPtr<LocalFromWorldPoint> LocalPointFinder;
00568 FieldSelector FieldSelection;
00569
00570 double step_size;
00571
00572 double scale;
00573
00574 bool valid_point;
00575
00576 pair<Eagle::PhysicalSpace::point, std::string>localPoint;
00577
00578 T CurrentMetric;
00579 tvector DataDir;
00580 RefPtr<ChristoffelXYZ<T> > Gamma;
00581
00582 enum { N = T::Dims * 2 };
00583 typedef double real;
00584
00585 GeodesicDifferentialEquation()
00586 {}
00587
00588 virtual ~GeodesicDifferentialEquation() {}
00589
00590 struct RealArray_t : FixedArray<real, N>
00591 {
00592 void init(int) {}
00593 };
00594
00595 struct DOPVarsArray_t : FixedArray<Traum::dop_vars<real>, N >
00596 {
00597 void init(int) {}
00598 };
00599
00600 void setData(const tvector&DataDirP, const RefPtr<LocalFromWorldPoint>&LocalPointFinderP, const FieldSelector&FieldSelectionP,
00601 const double step_sizeP, const double scaleP)
00602 {
00603 DataDir = DataDirP;
00604 LocalPointFinder = LocalPointFinderP;
00605 FieldSelection = FieldSelectionP;
00606 step_size = step_sizeP;
00607 scale = scaleP;
00608 }
00609
00610 void Accel(real s, const real *x, const real *v, real *d2x_ds2)
00611 {
00612 point Q;
00613 tvector q_dot;
00614
00615 for(int k=0; k < point::Dims; k++)
00616 {
00617 Q[k] = x[k];
00618 q_dot[k] = v[k];
00619 DataDir[k] = v[k];
00620 }
00621
00622
00623
00624 bool test = FieldLines::getChristoffelXYZ<T>( Q, LocalPointFinder, FieldSelection, step_size,
00625 Gamma,
00626 CurrentMetric,
00627 localPoint );
00628 if(!test)
00629 {
00630 puts("GeodesicDifferentialEquation::DifferentialEquation::Accel() No valid point found!");fflush(stdout);
00631 return;
00632 }
00633
00634 tvector qdd = getQddot<T>( Gamma, q_dot);
00635
00636 for(int k=0; k < point::Dims; k++)
00637 d2x_ds2[k] = qdd[k];
00638 }
00639
00640 };
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657 }
00658 #endif // __GEODESCICS_HPP_