FiberVISH 0.2
Fish - The Fiber Bundle API for the Vish Visualization Shell
Geodesics.hpp
1#ifndef __GEODESCICS_HPP_
2#define __GEODESCICS_HPP_
3
4#include <fish/fiber/vector/CubicIpol.hpp>
5#include <fish/fiber/vector/LinearIpol.hpp>
6#include <eagle/QuadraticMatrix.hpp>
7#include <eagle/Coordinates.hpp>
8#include <fish/fiber/baseop/LocalFromWorldPoint.hpp>
9#include <ode/dop853.hpp>
10
11//#define VERBOSE
12
13namespace Fiber
14{
15
16template <>
17struct LinearIpolZeroDerivativeTrait<metric33>
18{
19static metric33 zero()
20 {
21 metric33 m;
22 m.set(0.0);
23 return m;
24 }
25};
26
27template <>
28struct LinearIpolZeroDerivativeTrait<metric44>
29{
30static metric44 zero()
31 {
32 metric44 m;
33 m.set(0.0);
34 return m;
35 }
36};
37
38//template <>
39// struct CubicIpolZeroDerivativeTrait<metric33>
40// {
41// static metric33 zero()
42// {
43// metric33 m;
44// m.set(0.0);
45// return m;
46// }
47// };
48
49// template <>
50// struct CubicIpolZeroDerivativeTrait<metric44>
51// {
52// static metric44 zero()
53// {
54// metric44 m;
55// m.set(0.0);
56// return m;
57// }
58// };
59
60}
61
62
63namespace FieldLines
64{
65
66struct Geodesic{};
67
68struct GeodesicSpatial{};
69
70template<class T>
71double getDet(const Eagle::Quadratic<T::Dims,double>&q)
72{
73 return Eagle::Determinantor<T::Dims>::compute(q);
74}
75
76template<class T>
77class ChristoffelXYZ : public MemCore::ReferenceBase< ChristoffelXYZ<T> >
78{
79 typedef typename T::Chart_t Chart_t;
80 typedef typename Eagle::Coordinates<Chart_t>::point point;
81 // store metric tensor
82 T g;
83
84 Eagle::Quadratic<T::Dims, double> g_inv;
85
86 // store precomputed finite differences of g // array size 3
87 T g_diff[T::Dims];
88
89// double detG;
90
91 bool ok, normalizeInversion, doInversion;
92
93 // consts for access
94// enum { MU = 0, NU, LAMBDA };
95// enum { X = 0, Y, Z};
96
97public:
98 ChristoffelXYZ( const T&gP, const T&g_xP, const T&g_yP, const T&g_zP,
99 const double dxP, const double dyP, const double dzP,
100 bool normalizeInversionP, bool doInversionP )
101 : MemCore::ReferenceBase< ChristoffelXYZ >(this)
102 , g( gP )
103 , ok( true )
104 , normalizeInversion( normalizeInversionP )
105 , doInversion( doInversionP )
106 {
107 std::cout << "aber jetzt nit im ersnst!!!";
108 if(T::Dims == 3)
109 {
110 enum { X = 0, Y, Z};
111
112 if(dxP != 0.0)
113 for(unsigned i = 0; i < T::SIZE; i++)
114 g_diff[X][i] = ( g_xP[i] - g[i] ) / dxP;
115
116 if(dyP != 0.0)
117 for(unsigned i = 0; i < T::SIZE; i++)
118 g_diff[Y][i] = ( g_yP[i] - g[i] ) / dyP;
119
120 if(dzP != 0.0)
121 for(unsigned i = 0; i < T::SIZE; i++)
122 g_diff[Z][i] = ( g_zP[i] - g[i] ) / dzP;
123
124 }
125
126 if(T::Dims == 4)
127 {
128 enum { t = Eagle::STA::CartesianChart4D<>::T,
129 X = Eagle::STA::CartesianChart4D<>::X,
130 Y = Eagle::STA::CartesianChart4D<>::Y,
131 Z = Eagle::STA::CartesianChart4D<>::Z };
132
133 if(dxP != 0.0)
134 for(unsigned i = 0; i < T::SIZE; i++)
135 g_diff[X][i] = ( g_xP[i] - g[i] ) / dxP;
136
137 if(dyP != 0.0)
138 for(unsigned i = 0; i < T::SIZE; i++)
139 g_diff[Y][i] = ( g_yP[i] - g[i] ) / dyP;
140
141 if(dzP != 0.0)
142 for(unsigned i = 0; i < T::SIZE; i++)
143 g_diff[Z][i] = ( g_zP[i] - g[i] ) / dzP;
144
145 }
146 }
147
148 ChristoffelXYZ( const T&gP, const T&g_dxP, const T&g_dyP, const T&g_dzP,
149 bool normalizeInversionP, bool doInversionP )
150 : MemCore::ReferenceBase< ChristoffelXYZ<T> >(this)
151 , g(gP)
152 , ok(true)
153 , normalizeInversion( normalizeInversionP )
154 , doInversion( doInversionP )
155 {
156 //std::cout << g_dxP << endl << g_dyP << endl << g_dzP << endl;
157 //std::cout << "------";
158 if(T::Dims == 3)
159 {
160 enum { X = 0, Y, Z};
161 g_diff[X] = g_dxP;
162 g_diff[Y] = g_dyP;
163 g_diff[Z] = g_dzP;
164 }
165 if(T::Dims == 4)
166 {
167 enum { t = Eagle::STA::CartesianChart4D<>::T,
168 X = Eagle::STA::CartesianChart4D<>::X,
169 Y = Eagle::STA::CartesianChart4D<>::Y,
170 Z = Eagle::STA::CartesianChart4D<>::Z };
171
172 g_diff[t].set(0.0);
173 g_diff[X] = g_dxP;
174 g_diff[Y] = g_dyP;
175 g_diff[Z] = g_dzP;
176 }
177 try
178 {
179 Eagle::Quadratic<T::Dims, double> q (g);
180
181 double fac = q[0][0];
182
183 if( normalizeInversion )
184 {
185 for( size_t i = 0; i < T::Dims; i++ )
186 for( size_t j = 0; j < T::Dims; j++ )
187 q[i][j] = q[i][j] / fac;
188 }
189// std::cout << "q scaled: " << q << std::endl;
190
191 Eagle::ComputeInverse(g_inv, q);
192
193// std::cout << "g_inv: " << g_inv << std::endl;
194
195// Eagle::Matrix<T::Dims,T::Dims, double> &mg = q ;
196// Eagle::Matrix<T::Dims,T::Dims, double> &mg_i = g_inv;
197// Eagle::Matrix<T::Dims,T::Dims, double> I = mg * mg_i;
198
199// std::cout << "controll inversion: " << I << std::endl;
200 if( normalizeInversion )
201 {
202 for( size_t i = 0; i < T::Dims; i++ )
203 for( size_t j = 0; j < T::Dims; j++ )
204 g_inv[i][j] = g_inv[i][j] / fac;
205 }
206// std::cout << "g_inv scaled: " << g_inv << std::endl;
207
208 }
209 catch(const Eagle::DegeneratedMatrix&)
210 {
211 puts("ChristoffelXYZ<>::ChristoffelXYZ ERROR Degenerated Matrix in inversion."); fflush(stdout);
212 ok = false;
213 }
214 }
215
216 ChristoffelXYZ( const T&gP, const T&g_dxP, const T&g_dyP, const T&g_dzP, const T&g_dtP,
217 bool normalizeInversionP, bool doInversionP )
218 : MemCore::ReferenceBase< ChristoffelXYZ<T> >(this)
219 , g(gP)
220 , ok(true)
221 , normalizeInversion( normalizeInversionP )
222 , doInversion( doInversionP )
223 {
224 //std::cout << g_dxP << endl << g_dyP << endl << g_dzP << endl;
225 //std::cout << "------";
226 if(T::Dims == 3)
227 {
228 enum { X = 0, Y, Z};
229 g_diff[X] = g_dxP;
230 g_diff[Y] = g_dyP;
231 g_diff[Z] = g_dzP;
232 }
233 if(T::Dims == 4)
234 {
235 enum { t = Eagle::STA::CartesianChart4D<>::T,
236 X = Eagle::STA::CartesianChart4D<>::X,
237 Y = Eagle::STA::CartesianChart4D<>::Y,
238 Z = Eagle::STA::CartesianChart4D<>::Z };
239
240 g_diff[t] = g_dtP;
241 g_diff[X] = g_dxP;
242 g_diff[Y] = g_dyP;
243 g_diff[Z] = g_dzP;
244 }
245 try
246 {
247 Eagle::Quadratic<T::Dims, double> q (g);
248
249 double fac = q[0][0];
250
251 if( normalizeInversion )
252 {
253 for( size_t i = 0; i < T::Dims; i++ )
254 for( size_t j = 0; j < T::Dims; j++ )
255 q[i][j] = q[i][j] / fac;
256 }
257// std::cout << "q scaled: " << q << std::endl;
258
259 Eagle::ComputeInverse(g_inv, q);
260
261// std::cout << "g_inv: " << g_inv << std::endl;
262
263// Eagle::Matrix<T::Dims,T::Dims, double> &mg = q ;
264// Eagle::Matrix<T::Dims,T::Dims, double> &mg_i = g_inv;
265// Eagle::Matrix<T::Dims,T::Dims, double> I = mg * mg_i;
266
267// std::cout << "controll inversion: " << I << std::endl;
268 if( normalizeInversion )
269 {
270 for( size_t i = 0; i < T::Dims; i++ )
271 for( size_t j = 0; j < T::Dims; j++ )
272 g_inv[i][j] = g_inv[i][j] / fac;
273 }
274// std::cout << "g_inv scaled: " << g_inv << std::endl;
275
276 }
277 catch(const Eagle::DegeneratedMatrix&)
278 {
279 puts("ChristoffelXYZ<>::ChristoffelXYZ ERROR Degenerated Matrix in inversion."); fflush(stdout);
280 ok = false;
281 }
282 }
283
284 bool isOk(){ return ok;}
285
286 void getMetricfromQuadratic(const Eagle::Quadratic<T::Dims,double>&q, T&g)
287 {
288 for(unsigned i = 0; T::SIZE; i++)
289 g[i] = q[i];
290 }
291
292
293
294 // christoffel symbols function for on demand symbol computing, later there may be the need to
295 // store them and reuse them when possible.
296 const double operator()(const unsigned mu, const unsigned nu, const unsigned lambda) const
297 {
298 double gamma = 0.0;
299
300 if(mu > T::Dims-1 || nu > T::Dims-1 || lambda > T::Dims-1)
301 return 0.0;
302
303 if( doInversion )
304 for(unsigned i = 0; i < T::Dims; i++)
305 gamma += g_inv(lambda, i) * ( g_diff[nu](mu,i) + g_diff[mu](nu,i) - g_diff[i](mu, nu) );
306 else
307 for(unsigned i = 0; i < T::Dims; i++)
308 gamma += g(lambda, i) * ( g_diff[nu](mu,i) + g_diff[mu](nu,i) - g_diff[i](mu, nu) );
309
310 gamma /= 2;
311
312 return gamma;
313 }
314
315 const double gammaDebug(const unsigned mu, const unsigned nu, const unsigned lambda) const
316 {
317 double gamma = 0.0;
318
319 if(mu > T::Dims-1 || nu > T::Dims-1 || lambda > T::Dims-1)
320 return 0.0;
321 printf( "mu nu lam %d %d %d", mu, nu, lambda);
322 for(unsigned i = 0; i < T::Dims; i++)
323 {
324
325 double tmp = g_inv(lambda, i) * ( g_diff[nu](mu,i) + g_diff[mu](nu,i) - g_diff[i](mu, nu) );
326 gamma += tmp;
327 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);
328 }
329 gamma /= 2;
330 printf( "gamma: %f\n", gamma);
331 return gamma;
332 }
333
334
335 const T& getG() { return g; }
336
337 const T& getGDiff( unsigned iP ) { return g_diff[iP]; }
338
339 void speak() const
340 {
341
342 std::cout << "Christoffel content:\n" << "g " << endl << g << "g_inv: " << endl << g_inv;
343 std::cout << "normalize: " << normalizeInversion << std::endl;
344 std::cout << "invert: " << doInversion << std::endl;
345
346 if(T::Dims == 3)
347 std::cout << "g,x" << std::endl <<g_diff[0]
348 << "g,y" << std::endl <<g_diff[1]
349 << "g,z" << std::endl <<g_diff[2];
350
351 if(T::Dims == 4)
352 std::cout << "g,t" << std::endl <<g_diff[0]
353 << "g,x" << std::endl <<g_diff[1]
354 << "g,y" << std::endl <<g_diff[2]
355 << "g,z" << std::endl <<g_diff[3];
356
357 cout << "Gamma (T)|X|Y|Z|:" << endl;
358
359 for( unsigned lam = 0; lam < T::Dims; lam++)
360 {
361 for( unsigned mu = 0; mu < T::Dims; mu++)
362 {
363 for( unsigned nu = 0; nu < T::Dims; nu++)
364 {
365 printf("%f ", (*this)(mu,nu,lam));
366 }
367 printf("\n");
368 }
369 printf("-----------------------------\n");
370 }
371
372
373 // for( unsigned lam = 0; lam < T::Dims; lam++)
374 // {
375 // for( unsigned mu = 0; mu < T::Dims; mu++)
376 // {
377 // for( unsigned nu = 0; nu < T::Dims; nu++)
378 // {
379 // gammaDebug(mu, nu, lam);
380 // }
381 // printf("\n");
382 // }
383 // printf("-----------------------------\n");
384 // }
385 }
386
387};
388
389
392template<class T>
393typename Eagle::Coordinates<typename T::Chart_t >::vector getQddot( RefPtr< ChristoffelXYZ<T> >& Gamma,
394 typename Eagle::Coordinates<typename T::Chart_t >::vector&q_dot )
395{
396typename Eagle::Coordinates<typename T::Chart_t >::vector q_ddot;
397
398 q_ddot.set(0.0);
399
400 // q-ddot^lamba = - Gamma^lambda_mu_nu * q-dot^mu * q-dot^nu
401
402 //cout << "getQddot qdot:" << q_dot << endl;
403
404 for(index_t lambda = 0; lambda < T::Dims; lambda++)
405 for( index_t mu = 0; mu < T::Dims; mu++)
406 for(index_t nu = 0; nu < T::Dims; nu++)
407 q_ddot[lambda] -= (*Gamma)( mu, nu, lambda ) * q_dot[mu] * q_dot[nu];
408
409 //cout << "getQddot qddot:" << q_ddot << endl;
410
411 return q_ddot;
412}
413
414
415template<class TT>
416typename Eagle::Coordinates<typename TT::Chart_t >::vector getQdot( const typename Eagle::Coordinates<typename TT::Chart_t >::vector& q_dot,
417 const TT&g, const bool big = true)
418{
419 if (TT::Dims==3)
420 return q_dot;
421
422 enum { T = STA::CartesianChart4D<>::T,
423 X = STA::CartesianChart4D<>::X,
424 Y = STA::CartesianChart4D<>::Y,
425 Z = STA::CartesianChart4D<>::Z };
426
427const double eps = 1.0e-4;
428
429 //std::cout << "qdot before:" << q_dot << std::endl;
430
431typename Eagle::Coordinates<typename TT::Chart_t >::vector v = q_dot;
432
433 if( !isnormal(v[X]) ) v[X] = 0.0;
434 if( !isnormal(v[Y]) ) v[Y] = 0.0;
435 if( !isnormal(v[Z]) ) v[Z] = 0.0;
436
437
438 if(v[X] < eps && v[X] > -eps &&
439 v[Y] < eps && v[Y] > -eps &&
440 v[Z] < eps && v[Z] > -eps)
441 {
442 puts("Creating Initial Condition for falling geodesic.");
443 v[T] = 1.0;
444 return v;
445 }
446
447 if(g(T,T) == 0.0)
448 {
449 puts("getQdot::Error g(T,T) == 0! Division by 0.");
450 return v;
451 }
452
453 //puts("hm");
454
455double A = 2 * ( g(X,T)*v[X] + g(Y,T)*v[Y] + g(Z,T)*v[Z] );
456double B = g(X,X)*v[X]*v[X] + g(Y,Y)*v[Y]*v[Y] + g(Z,Z)*v[Z]*v[Z]+
457 2*( g(X,Y)*v[X]*v[Y] + g(X,Z)*v[X]*v[Z] + g(Y,Z)*v[Y]*v[Z] );
458
459double p = A/g(T,T);
460double q= B/g(T,T);
461
462double w = p*p/4 - q;
463
464 if(w < 0.0)
465 {
466 puts("getQdot::Error direction initialization srqt < 0!!!!");
467 return v;
468 }
469
470 double sw = sqrt( w );
471
472 double v1 = -p/2 + sw;
473 double v2 = -p/2 - sw;
474
475 if(big)
476 v[T] = (v1 > v2)? v1 : v2; // use bigger value of the solution of the quadratic equation
477 else
478 v[T] = (v1 > v2)? v2 : v1;
479
480 return v;
481}
482
483
486template<class T>
487bool getChristoffelXYZ( const typename Eagle::Coordinates<typename T::Chart_t >::point& CurrentPoint,
488 const MemCore::RefPtr<LocalFromWorldPoint>&LocalPointFinder,
489 const FieldSelector&FieldSelection, const double step_size,
490 RefPtr< ChristoffelXYZ<T> >& Christoffels,
491 T& CurrentMetric,
492 LocalPoint& localPoint,
493 bool normalizeInversion, bool doInversion )
494 {
495#ifdef VERBOSE
496 puts("getChristoffelXYZ<class T>() enter" );
497 cout << "AtomicIntegrator::doEuler() invNorm: " << normalizeInversion << endl;
498 cout << "AtomicIntegrator::doEuler() invEnable: " << doInversion << endl;
499#endif
500 typedef typename T::Chart_t Chart_t;
501 typedef typename Coordinates<Chart_t>::point point;
502
503 // get spatial point
504 Eagle::PhysicalSpace::point current (CurrentPoint[0+T::Dims-3], CurrentPoint[1+T::Dims-3], CurrentPoint[2+T::Dims-3]);
505
506
507#ifdef VERBOSE
508 puts("getChristoffelXYZ<class T>() getting local point" );
509#endif
510 bool test = LocalPointFinder->get( current, localPoint );
511
512 if(!test)
513 {
514 puts("getChristoffelXYZ()1 error finding local point at:");
515 cout << current << endl;
516 return false;
517 }
518
519 RefPtr<FragmentID> FragName = localPoint.frag_id;
520 const Eagle::PhysicalSpace::point &float_index = localPoint.local_point;
521
522
523#ifdef VERBOSE
524 std::cout << "getChristoffelXYZ<class T>() getting field: " << FieldSelection.getFieldName() << std::endl;
525#endif
526
527 RefPtr<Field> Field;
528 Field = FieldSelection.getField();
529
530#ifdef VERBOSE
531 cout << "FieldLines::getChristoffelXYZ local point: " << float_index;
532 cout << "in Fragment:" << FragName->Name() << endl;
533 cout << "getting Field: " << FieldSelection.getFieldName() << endl;
534#endif
535 if(!Field)
536 {
537 puts("AtomicIntegrator<metric33, Geodesic>doit2(); ERROR: No Field0 found in Time iteration");
538 return false;
539 }
540
541
542 RefPtr<FragmentID> FragID = FragName;
543// if (FragName.compare("") != 0 )
544// FragID = new FragmentID( FragName );
545
546 RefPtr<CreativeArrayBase> cBase = Field->getCreator( FragID );
547 assert(cBase);
548
549 RefPtr<MemBase> mBase = cBase->create();
550 assert(mBase);
551
552
553
554 RefPtr<MemArray<3, T> > ToIntArr0 = mBase;
555 ToIntArr0 = Field->getCreator( FragID )->create();
556
557 if(!ToIntArr0)
558 {
559 cout << "getChristoffelXYZ<T> dims:" << T::Dims << endl;
560 mBase.speak("getChristoffelXYZ<T> ------>");
561
562 puts("getChristoffelXYZ()<T>:ERROR retrieving matching field. Maybe Dims of instantiated metric and Dims of field data does not match. 33 44?");
563 return false;
564 }
565
566 assert( ToIntArr0 );
567
568
569
570 // NOTE: Here optionally use another interpolator, eventually cubic, or
571 // self-aligning for eigenvector fields.
572 //
573
574// cout << "before interpol" << endl;
575
576 T g_x, g_y, g_z;
577 T g_xl, g_xp, g_xxp, g_yp, g_zp, g_yzp, g_yyzp;
578
579 Interpolate<3, T, CubicIpol<T> > myField0 ( *ToIntArr0, float_index);
580 Interpolate<3, T, CubicIpol<T>, double, NoDelimiter<T>, 0 > myField0_x( *ToIntArr0, float_index);
581 Interpolate<3, T, CubicIpol<T>, double, NoDelimiter<T>, 1 > myField0_y( *ToIntArr0, float_index);
582 Interpolate<3, T, CubicIpol<T>, double, NoDelimiter<T>, 2 > myField0_z( *ToIntArr0, float_index);
583 // Interpolate<3, T, LinearIpol<T> > myField0 ( *ToIntArr0, float_index);
584// Interpolate<3, T, LinearIpol<T>, double, NoDelimiter<T>, 0 > myField0_xl( *ToIntArr0, float_index);
585 // Interpolate<3, T, LinearIpol<T>, double, NoDelimiter<T>, 1 > myField0_y( *ToIntArr0, float_index);
586 // Interpolate<3, T, LinearIpol<T>, double, NoDelimiter<T>, 2 > myField0_z( *ToIntArr0, float_index);
587
588 CurrentMetric = myField0.eval();
589 g_x = myField0_x.eval();
590 g_y = myField0_y.eval();
591 g_z = myField0_z.eval();
592
593// g_xl= myField0_xl.eval();
594
595 // double eps = -1;
596
597 // Eagle::point3 xp = float_index + Eagle::tvector3(eps,0,0 );
598 // Eagle::point3 xxp = float_index + Eagle::tvector3(-eps,0,0 );
599
600 // Eagle::point3 yp = float_index + Eagle::tvector3(0,eps,0 );
601 // Eagle::point3 zp = float_index + Eagle::tvector3(0,0,eps );
602 // Eagle::point3 yzp = float_index + Eagle::tvector3(0,0,eps ) + Eagle::tvector3(0,eps,0 );
603 // Eagle::point3 yyzp = float_index + Eagle::tvector3(0,0,eps ) - Eagle::tvector3(0,eps,0 );
604
605 //Interpolate<3, T, LinearIpol<T> > myField0_xp( *ToIntArr0, xp );
606 //Interpolate<3, T, LinearIpol<T> > myField0_xxp( *ToIntArr0, xxp );
607 //Interpolate<3, T, LinearIpol<T> > myField0_yp( *ToIntArr0, yp );
608 //Interpolate<3, T, LinearIpol<T> > myField0_zp( *ToIntArr0, zp );
609 // Interpolate<3, T, LinearIpol<T> > myField0_yzp( *ToIntArr0, yzp );
610 // Interpolate<3, T, LinearIpol<T> > myField0_yyzp( *ToIntArr0, yyzp );
611
612 // g_xp = myField0_xp.eval();
613 // g_xxp = myField0_xxp.eval();
614
615 // g_yp = myField0_yp.eval();
616 // g_zp = myField0_zp.eval();
617 // g_yzp = myField0_yzp.eval();
618 // g_yyzp = myField0_yyzp.eval();
619
620 RefPtr<Fiber::Field> field = FieldSelection.getCartesianPositions();
621 RefPtr<Fiber::LinearDirectProductMemArray<Eagle::point3> > FData = field->getData();
622 Eagle::point3 delta = FData->Delta();
623
624
625 Eagle::Quadratic<T::Dims, double> q (CurrentMetric);
626 double det = getDet<T>(q);
627
628
629
630 if(det < 0.0 && T::Dims == 3) // check if det gets negative
631 {
632 cout << "g:" << q << endl;
633 cout << "det: " << det << endl;
634 puts("getChristoffelXYZ()<T>:warning, Det of G < 0.");
635 return false;
636 }
637
638 g_x /= delta[0]; // normalize index derivative by interval length
639 g_xl /= delta[0];
640 g_y /= delta[1];
641 g_z /= delta[2];
642
643#ifdef VERBOSE
644/*
645 cout << "det(g): " << det << endl;
646 cout << "delta: " << delta << endl;
647 cout << "aus dem interpolator:" << endl << CurrentMetric << g_x << g_y << g_z;
648
649 cout << endl <<"g:" << endl
650 << float_index << endl << CurrentMetric << "variiert xyz:" << endl
651 << xp << endl << g_xp
652 << xxp << endl << g_xxp
653
654 << yp << endl << g_yp
655 << zp << endl << g_zp
656
657 << "g,x lin" << endl << g_xl;
658 // << yzp << endl << g_yzp
659 // << yyzp << endl << g_yyzp ;
660*/
661#endif
662 Christoffels = new ChristoffelXYZ<T>( CurrentMetric, g_x, g_y, g_z, normalizeInversion, doInversion );
663
664#ifdef VERBOSE
665 std::cout << "Hilfe Christoffel wo bist du?" << std::endl;
666 cout << "current point: " << CurrentPoint;
667 Christoffels->speak();
668#endif
669 if( Christoffels->isOk() == false)
670 {
671 puts("getChristoffelXYZ()<T>:warning, christoffels probably degenerated matrix.");
672 std::cout << "Hilfe Christoffel wo bist du?" << std::endl;
673 cout << "current point: " << CurrentPoint;
674 Christoffels->speak();
675 return false;
676 }
677
678 return true;
679 }
680
681
682
684// different function signature. A bit cleaned and used in front_geodesics, still lots of work can be done to improve ...
685// todo: curvilinear coordinates
686// distinction of type of coordinate could be done somewhere else? not for every christoffel call?
688template<class T>
689bool getChristoffelXYZ( const typename Eagle::Coordinates<typename T::Chart_t >::point& current_pointP,
690 const LocalPoint& local_pointP,
691 const LocalFromWorldPoint finderP,
692 const RefPtr<Field>& tensor_fieldP,
693 const RefPtr<Field>& coord_fieldP,
694 RefPtr< ChristoffelXYZ<T> >& christoffelsP,
695 T& current_metricP,
696 bool normalizeInversion, bool doInversion, bool central_derivative = false )
697 {
698#ifdef VERBOSE
699 puts("getChristoffelXYZ<class T>() enter" );
700 cout << "AtomicIntegrator::doEuler() invNorm: " << normalizeInversion << endl;
701 cout << "AtomicIntegrator::doEuler() invEnable: " << doInversion << endl;
702#endif
703 typedef typename T::Chart_t Chart_t;
704 typedef typename Coordinates<Chart_t>::point point;
705
706 Eagle::Quadratic<T::Dims, double>& q (current_metricP);
707 double det = getDet<T>(q);
708
709 if(det < 0.0 && T::Dims == 3) // check if det gets negative
710 {
711 cout << "g:" << q << endl;
712 cout << "det: " << det << endl;
713 puts("getChristoffelXYZ()<T>:warning, Det of G < 0.");
714 return false;
715 }
716
717 Eagle::PhysicalSpace::point current ( current_pointP[0+T::Dims-3],
718 current_pointP[1+T::Dims-3],
719 current_pointP[2+T::Dims-3] );
720#ifdef VERBOSE
721 puts("getChristoffelXYZ<class T>() getting local point" );
722#endif
723
724 RefPtr<FragmentID> frag_name = local_pointP.frag_id;
725 const Eagle::PhysicalSpace::point &float_index = local_pointP.local_point;
726
727#ifdef VERBOSE
728 cout << "FieldLines::getChristoffelXYZ local point: " << float_index;
729 cout << "in Fragment:" << frag_name->Name() << endl;
730#endif
731 if(!tensor_fieldP)
732 {
733 puts("AtomicIntegrator<metric33, Geodesic>doit2(); ERROR: No Field0 found in Time iteration");
734 return false;
735 }
736
737 RefPtr<FragmentID> frag_ID = frag_name;
738
739 RefPtr<MemArray<3, T> > ToIntArr0 = tensor_fieldP->getData( frag_ID );
740
741 // having uniform or curvilinear data or rectilinear data
742 if( ToIntArr0 )
743 {
744 T g_x, g_y, g_z;
745 T g_xl, g_xp, g_xxp, g_yp, g_zp, g_yzp, g_yyzp;
746
747 Interpolate<3, T, CubicIpol<T> > myField0 ( *ToIntArr0, float_index);
748 current_metricP = myField0.eval();
749
750 RefPtr<Fiber::LinearDirectProductMemArray<Eagle::point3> > FData = coord_fieldP->getData();
751
752 // having uniform grid
753 if( FData )
754 {
755 Interpolate<3, T, CubicIpol<T>, double, NoDelimiter<T>, 0 > myField0_x( *ToIntArr0, float_index);
756 Interpolate<3, T, CubicIpol<T>, double, NoDelimiter<T>, 1 > myField0_y( *ToIntArr0, float_index);
757 Interpolate<3, T, CubicIpol<T>, double, NoDelimiter<T>, 2 > myField0_z( *ToIntArr0, float_index);
758 g_x = myField0_x.eval();
759 g_y = myField0_y.eval();
760 g_z = myField0_z.eval();
761
762 Eagle::point3 delta = FData->Delta();
763
764 Eagle::Quadratic<T::Dims, double> q (current_metricP);
765 double det = getDet<T>(q);
766
767 if(det < 0.0 && T::Dims == 3) // check if det gets negative
768 {
769 cout << "g:" << q << endl;
770 cout << "det: " << det << endl;
771 puts("getChristoffelXYZ()<T>:warning, Det of G < 0.");
772 return false;
773 }
774
775 g_x /= delta[0]; // normalize index derivative by interval length
776 g_xl /= delta[0];
777 g_y /= delta[1];
778 g_z /= delta[2];
779
780#ifdef VERBOSE
781 cout << "det(g): " << det << endl;
782 cout << "delta: " << delta << endl;
783 cout << "aus dem interpolator:" << endl << CurrentMetric << g_x << g_y << g_z;
784
785 cout << endl <<"g:" << endl
786 << float_index << endl << CurrentMetric << "variiert xyz:" << endl
787 << xp << endl << g_xp
788 << xxp << endl << g_xxp
789
790 << yp << endl << g_yp
791 << zp << endl << g_zp
792
793 << "g,x lin" << endl << g_xl;
794 // << yzp << endl << g_yzp
795 // << yyzp << endl << g_yyzp ;
796#endif
797 christoffelsP = new ChristoffelXYZ<T>( current_metricP, g_x, g_y, g_z, normalizeInversion, doInversion );
798
799
800#ifdef VERBOSE
801 std::cout << "Hilfe Christoffel wo bist du?" << std::endl;
802// cout << "current point: " << CurrentPoint;
803 christoffelsP->speak();
804#endif
805 if( christoffelsP->isOk() == false)
806 {
807 puts("getChristoffelXYZ()<T>:warning, christoffels probably degenerated matrix.");
808 std::cout << "Hilfe Christoffel wo bist du?" << std::endl;
809 cout << "current point: " << current_pointP;
810 christoffelsP->speak();
811 return false;
812 }
813
814 return true;
815 }
816
817 // curvilinear data
818 RefPtr<MemArray<3, Eagle::point3> > f_data_curvi = coord_fieldP->getData();
819 if( f_data_curvi )
820 {
821 assert( local_pointP.cell_size.x() > 0.0 );
822
823 // do a finite difference (x1 - x0) / h in world coordinates with h the cell size
824
825 LocalPoint local_x( local_pointP.world_point + tvector( local_pointP.cell_size.x(), 0, 0 ) );
826 LocalPoint local_y( local_pointP.world_point + tvector( 0, local_pointP.cell_size.y(), 0 ) );
827 LocalPoint local_z( local_pointP.world_point + tvector( 0, 0, local_pointP.cell_size.z() ) );
828 if( !finderP.getViaTree( local_x ) || !finderP.getViaTree( local_y ) || !finderP.getViaTree( local_z ) )
829 {
830 puts("getChristoffelXYZ()<T>:warning, local point not found for x+h");
831 return false;
832 }
833
834
835
836 Interpolate<3, T, CubicIpol<T> > myField0x( *ToIntArr0, local_x );
837 Interpolate<3, T, CubicIpol<T> > myField0y( *ToIntArr0, local_y );
838 Interpolate<3, T, CubicIpol<T> > myField0z( *ToIntArr0, local_z );
839
840 g_x = myField0x.eval();
841 g_y = myField0y.eval();
842 g_z = myField0z.eval();
843
844
845
846 if( central_derivative == true)
847 {
848 LocalPoint local_x_( local_pointP.world_point + tvector( -local_pointP.cell_size.x(), 0, 0 ) );
849 LocalPoint local_y_( local_pointP.world_point + tvector( 0,-local_pointP.cell_size.y(), 0 ) );
850 LocalPoint local_z_( local_pointP.world_point + tvector( 0,0,-local_pointP.cell_size.z() ) );
851 if( !finderP.getViaTree( local_x_ ) || !finderP.getViaTree( local_y_ ) || !finderP.getViaTree( local_z_ ))
852 {
853 puts("getChristoffelXYZ()<T>:warning, local point not found for x-h");
854 return false;
855 }
856 Interpolate<3, T, CubicIpol<T> > myField0x_( *ToIntArr0, local_x_ );
857 Interpolate<3, T, CubicIpol<T> > myField0y_( *ToIntArr0, local_y_ );
858 Interpolate<3, T, CubicIpol<T> > myField0z_( *ToIntArr0, local_z_ );
859
860 //
861 g_x -= myField0x_.eval();
862 g_y -= myField0x_.eval();
863 g_z -= myField0x_.eval();
864
865 g_x /= (2*local_pointP.cell_size.x());
866 g_y /= (2*local_pointP.cell_size.y());
867 g_z /= (2*local_pointP.cell_size.z());
868 }
869 else
870 {
871 g_x = (g_x - current_metricP) / local_pointP.cell_size.x();
872 g_y = (g_y - current_metricP) / local_pointP.cell_size.y();
873 g_z = (g_z - current_metricP) / local_pointP.cell_size.z();
874 }
875
876 christoffelsP = new ChristoffelXYZ<T>( current_metricP, g_x, g_y, g_z, normalizeInversion, doInversion );
877
878#ifdef VERBOSE
879 std::cout << "Hilfe Christoffel wo bist du?" << std::endl;
880// cout << "current point: " << CurrentPoint;
881 christoffelsP->speak();
882#endif
883 if( christoffelsP->isOk() == false)
884 {
885 puts("getChristoffelXYZ()<T>:warning, christoffels probably degenerated matrix.");
886 std::cout << "Hilfe Christoffel wo bist du?" << std::endl;
887 cout << "current point: " << current_pointP;
888 christoffelsP->speak();
889 return false;
890 }
891
892 return true;
893
894 }
895
896 }
897
898 // fem or hexahedral data, or point cloud data
899 //RefPtr<MemArray<1, Eagle::point3> > f_data_curvi = coord_fieldP->getData();
900}
901
902
903
904template<class T>
905struct GeodesicDifferentialEquation : public Traum::VirtualSecondOrderDiffEquation <double>
906{
907typedef typename T::Chart_t Chart_t;
908typedef typename Eagle::Coordinates<Chart_t>::point point;
909typedef typename Eagle::Coordinates<Chart_t>::vector tvector;
910
911RefPtr<LocalFromWorldPoint> LocalPointFinder;
912FieldSelector FieldSelection;
913
914double step_size; // depricated, define stepsize for finite differential derivation of the metric tensors used for ChristoffelsXYZ
915
916
917double scale;
918
919bool valid_point;
920
921pair<Eagle::PhysicalSpace::point, RefPtr<FragmentID> > localPoint;
922
923T CurrentMetric;
924tvector DataDir;
925RefPtr<ChristoffelXYZ<T> > Gamma;
926
927bool normalizeInversion, doInversion;
928
929 enum { N = T::Dims * 2 };
930 typedef double real;
931
932 GeodesicDifferentialEquation()
933 {}
934
935 virtual ~GeodesicDifferentialEquation() {}
936
937 struct RealArray_t : FixedArray<real, N>
938 {
939 void init(int) {}
940 };
941
942 struct DOPVarsArray_t : FixedArray<Traum::dop_vars<real>, N >
943 {
944 void init(int) {}
945 };
946
947 void setData(const tvector&DataDirP, const RefPtr<LocalFromWorldPoint>&LocalPointFinderP, const FieldSelector&FieldSelectionP,
948 const double step_sizeP, const double scaleP)
949 {
950 DataDir = DataDirP;
951 LocalPointFinder = LocalPointFinderP;
952 FieldSelection = FieldSelectionP;
953 step_size = step_sizeP;
954 scale = scaleP;
955 }
956
957 void Accel(real s, const real *x, const real *v, real *d2x_ds2)
958 {
959 point Q;
960 tvector q_dot;
961
962 for(int k=0; k < point::Dims; k++)
963 {
964 Q[k] = x[k];
965 q_dot[k] = v[k];
966 DataDir[k] = v[k];
967 }
968// printf("Q: %f, %f, %f\n", Q[0], Q[1], Q[2] );
969// printf("V: %f, %f, %f\n", v[0], v[1], v[2] );
970
971 bool test = FieldLines::getChristoffelXYZ<T>( Q, LocalPointFinder, FieldSelection, step_size,
972 Gamma,
973 CurrentMetric,
974 localPoint, normalizeInversion, doInversion );
975 if(!test)
976 {
977 puts("GeodesicDifferentialEquation::DifferentialEquation::Accel() No valid point found!");fflush(stdout);
978 return;
979 }
980
981 tvector qdd = getQddot<T>( Gamma, q_dot);
982
983 for(int k=0; k < point::Dims; k++)
984 d2x_ds2[k] = qdd[k];
985 }
986
987};
988
989// double linearInterpolate(const double t0, const double d0, const double t1, const double d1, const double t)
990// {
991// return d0 + (d1 - d0) / (t1 - t0) * (t - t0);
992// }
993
994// metric33 linearInterpolate(const double t0, const metric33&m0, const double t1, const metric33&m1, const double t )
995// {
996// metric33 tmp;
997// for(unsigned i = 0; i < metric33::SIZE; i++)
998// tmp[i] = linearInterpolate(t0, m0[i], t1, m1[i], t);
999
1000// return tmp;
1001// }
1002
1003
1004}
1005#endif // __GEODESCICS_HPP_
sqrt
complex< _Tp > sqrt(const complex< _Tp > &)
endl
basic_ostream< _CharT, _Traits > & endl(basic_ostream< _CharT, _Traits > &__os)
std::cout
ostream cout
Fiber::CreativeIterator
An iterator with an optional DataCreator, which is just a class to intercept creation of data along a...
Definition CreativeIterator.hpp:34
Fiber::FieldSelection
An abstract selection of fields, that is given by names of fields and possible types for each field.
Definition FieldSelection.hpp:23
Fiber::FieldSelection::getFieldName
string getFieldName() const
Retrieve the name of a unique field, if only one is stored here.
Definition FieldSelection.cpp:68
Fiber::FieldSelector
An internal class that stores a couple of textual names.
Definition FieldSelector.hpp:18
Fiber::Field
A Field is a collection of CreativeArrayBase reference pointers which are accessed via FragmentID obj...
Definition Field.hpp:245
Fiber::Field::getCreator
RefPtr< CreativeArrayBase > getCreator(const RefPtr< FragmentID > &=nullptr) const
Get the creator of an unfragmented field.
Definition Field.cpp:603
Fiber::Iterator::set
void set(const T &Value)
Set all elements to the same value.
Definition vector/Iterator.hpp:724
Fiber::LocalFromWorldPoint
Use this class to compute a local cell index and fragment ID (if neccessary) of any coordinate field ...
Definition LocalFromWorldPoint.hpp:71
FieldLines::ChristoffelXYZ
Definition Geodesics.hpp:78
Eagle::Determinantor
Fiber
Given a fragmented field of curvilinear coordinates, (3D array of coordinates), build a uniform Grid ...
Definition FAQ.dox:2
Fiber::LinearIpolZeroDerivativeTrait
Helper class (type trait) for class LinerIpol, which for a given type T provides a static member func...
Definition LinearIpol.hpp:25
Fiber::LocalPoint
Definition LocalFromWorldPoint.hpp:39
FieldLines::GeodesicDifferentialEquation::DOPVarsArray_t
Definition Geodesics.hpp:943
FieldLines::GeodesicDifferentialEquation::RealArray_t
Definition Geodesics.hpp:938
FieldLines::GeodesicDifferentialEquation
Definition Geodesics.hpp:906
FieldLines::GeodesicSpatial
Definition Geodesics.hpp:68
FieldLines::Geodesic
Definition Geodesics.hpp:66
Generated on Sun Mar 1 2026 10:23:23 for FiberVISH by doxygen 1.9.8