IntegralLine.hpp

#include "IntegralFace.hpp"
#include <fish/lakeview/base/FieldInterpolator.hpp>

#include <ocean/plankton/VCreator.hpp>
#include <ocean/eagle/PhysicalSpace.hpp>
#include <ocean/eagle/STA.hpp>
#include <ocean/shrimp/VEnum.hpp>
#include <ocean/shrimp/PhysicalSpace.hpp>

#include <ocean/plankton/VTimer.hpp>
 
#include <bone/FishField.hpp>
#include <bone/FishSlice.hpp>

#include <ocean/shrimp/VObjectStatus.hpp>
#include <fish/fiber/baseop/LocalFromWorldPoint.hpp>
#include <fish/fiber/baseop/ExpandBBox.hpp>

#include <ode/dop853def.hpp>
#include <baseop/TangentialVectors.hpp>

#include <ocean/plankton/VTimer.hpp>

#include <grid/types/LineSet.hpp>




//#define VERBOSE

#define TIMEIT

namespace FieldLines
{

// using namespace Wizt;
// using namespace Fiber;
// using namespace Eagle;
// using namespace Eagle::PhysicalSpace;
// using namespace Traum;

typedef std::string FragmentID_t;

template<class FieldType>
struct  IntegrationPoint
{
        typedef typename FieldType::Chart_t Chart_t;
        typedef typename Coordinates<Chart_t>::point point;
        typedef typename Coordinates<Chart_t>::vector tvector;
        
        //point      location;
        point location;
        FieldType            data;
        Eagle::PhysicalSpace::point index_location;

        FragmentID_t fragment;

        bool is_end;

        tvector      direction;

        double length;

        IntegrationPoint()
        :is_end(false)
        {}
        void setLocation(point&locationP)
        {
                location = locationP;
        }

        void setDirection(tvector directionP)
        {
                direction = directionP;
        }

        void setData( const FieldType&dataP, const Eagle::PhysicalSpace::point&index_locationP, const FragmentID_t&fragmentP, const double lengthP )
        {
                data = dataP; 
                index_location = index_locationP; 
                fragment = fragmentP; 
                length = lengthP;
        }

        void setEnd()
        {
                is_end = true;
        }
};

// template<class P>
// struct Pointilizer
// {
//      typedef typename P::Chart_t point;
// };

// template<>
// class Pointilizer<Coordinates<Eagle::CartesianChart3D > >
// {
//      typedef Eagle::point3 point;
// };

// template<>
// class Pointilizer<Coordinates<STA::CartesianChart4D<> > >
// {
//      typedef point4 point;
// };

 /*
  Note: This kind of data structure is inefficient, it should 
  definitely be revised. With this layout, it requires copying data
  for getting these data into the fiber bundle layout.
 */
//typedef       RefPtr<Chunk<IntegrationPoint<metric33> > > IntegrationResult_t;

//      IntegrationResult_t x;


/*
This is the structure how it should be. This one is fiber-bundle compatible.

template<class FieldType>
struct  IntegrationResult
{
        RefPtr<Chunk<point> >        location;

        RefPtr<Chunk<FieldType> >            direction;

        RefPtr<Chunk<point> >        index_location;

        RefPtr<Chunk<FragmentID_t> > fragment;
};
*/


/*
                if( interpol == 0 )     
                {
                        VectorLinearInterpolator_t myField(*ToIntArr, float_index); 
                                InterpolData = myField.eval();
                }       
                else if( interpol == 1 )
                {
                        VectorCubicInterpolator_t myField(*ToIntArr, float_index); 
                                InterpolData = myField.eval();
                } 
                else
                {
                        if ( MemCore::InterfacePtr<Eqn_t> AnalyticField = *ToIntegrateField )
                        {
                                Eqn_t&AF = *AnalyticField; 
                                        InterpolData = AF(time, location);
                        } 
                        else
                        {
                        VectorLinearInterpolator_t myField(*ToIntArr, float_index); 
                                InterpolData = myField.eval();
                        }
                }
*/


struct AtomicDataBase
{
        RefPtr<LocalFromWorldPoint> LocalPointFinder;
        FieldSelector FieldSelection;
        RefPtr<BoundingBox> BB;
        double step_size; 
        double scale;
        bool inited;
//      unsigned interpol;

        AtomicDataBase( const RefPtr<LocalFromWorldPoint>&LocalPointFinderP, 
                        const FieldSelector&FieldSelectionP, const RefPtr<BoundingBox>&BBP, 
                        const double step_sizeP , const double scale_facP ) //, const unsigned interpolP )
        : LocalPointFinder(LocalPointFinderP)
        , FieldSelection(FieldSelectionP) // somehow this was not sufficient ?!? see below
        , BB(BBP)
        , step_size(step_sizeP)
        , scale(scale_facP)
        , inited(false)
          //, interpol(interpolP)
        {
                FieldSelection = FieldSelectionP;       
#ifdef VERBOSE
                
                cout << "AtomicDataBase::AtomicDataBase FieldSelectionP: " << FieldSelectionP.getFieldName() << endl; 
                cout << "AtomicDataBase::AtomicDataBase FieldSelection: " << FieldSelection.getFieldName() << endl; 
#endif          
        }
};


class IntegrationFrontContainerBase
{
protected:
        unsigned size;
        unsigned last_index;
        
        unsigned front_index;

public:
        IntegrationFrontContainerBase(const unsigned sizeP = 0)
        : size(sizeP)
        , last_index(0)
        , front_index(0)
                {}

        void setSize(const unsigned sizeP)     { size = sizeP; }
        unsigned getSize() const               { return size; }
        void setLastIndex(const unsigned last) { last_index = last; }
        void advanceIndex()                    { last_index = last_index++ % size; }
        unsigned getLastIndex() const          { return last_index; }  
        void advanceFrontIndex()               { front_index++; }
        unsigned getFrontIndex() const         { return front_index; }
};


template<typename FieldType, typename LineType>
class IntegrationFrontContainter
{

public:

        void insertDataToBundle( FieldSelector&FieldSelection, const string module_name )
        {
                puts("IntegrationFrontContainter<>::insertDataToBundle() to be implemented!");
        }
        
};




// type trait for atomic integration requires 1 Euler and 2 Dop853 functions
template<typename FieldType, typename LineType, int InterpolType>
class AtomicIntegrator
{
        typedef RefPtr<Chunk<IntegrationPoint<FieldType> > > IntegrationResult_t;
        typedef MemArray<3, FieldType> ToIntegrateArray_t; 
public:
        bool doEuler( std::vector<IntegrationPoint<FieldType> >&IntegrationLine, const double&time )
        {
                puts("FieldLines::AtomicIntegrator<FieldType>.doEuler() ERROR: Unsupported type FieldType doing nothing!");fflush(stdout); 
                return false;
        }

        // bool doEuler( IntegrationContainer<FieldType, LineType>&IntegrationData
        //            const double&time )

        bool doDop853( std::vector<IntegrationPoint<FieldType> >&IntegrationLine, 
                       const RefPtr<LocalFromWorldPoint>&LocalPointFinder,
                       const FieldSelector&FieldSelection,
                       const RefPtr<BoundingBox>&BB,
                       const double step_size, const double time,
                       const unsigned line_nr )
        {
                puts("FieldLines::AtomicIntegrator<FieldType>.doDop853() ERROR: Unsupported type FieldType doing nothing!");fflush(stdout); 
                return false;   
        }

        bool initDop853( const int number_of_lines )
        {
                puts("FieldLines::AtomicIntegrator<FieldType>.initDop853() ERROR: Unsupported type FieldType doing nothing!");fflush(stdout); 
                return false;
        }       

        void finalize( std::vector<IntegrationPoint<FieldType> >&IntegrationLine, const double&time )
        {
                puts("FieldLines::AtomicIntegrator<FieldType>.finalize() INFO: Unsupported type FieldType doing nothing!");fflush(stdout); 
                return;
        }
};



template <class Extractor, class Container>
RefPtr<MemBase> CollectLines(const Extractor&E, const std::vector<RefPtr<Chunk<Container> > >&Lines, size_t nPoints)
{
typedef typename Extractor::output_t output_t;
typedef MemArray<1, output_t> ResultMemArray_t;

size_t  nLines = Lines.size(); 
Ref<ResultMemArray_t> ResultArray( MIndex(nPoints) );
std::vector<output_t>&LineResult = ResultArray->myChunk()->get_vector(); 

index_t VertexNr = 0;
        for(index_t LineNr = 0; LineNr < nLines; LineNr++)
        {
        const std::vector<Container>&PointsPerLine = Lines[LineNr]->std_vector(); 
                for(index_t LineVertex = 0; LineVertex < PointsPerLine.size() ; LineVertex++ )
                {
                        LineResult[ VertexNr ] = E(PointsPerLine[ LineVertex ]); 
                        VertexNr++;
                }
        }
        return ResultArray;
}

template <class Container>
void CollectLineData(Representation&CartesianVertices, const std::vector<RefPtr<Chunk<Container> > >&Lines, size_t nPoints);


struct  CopyExtractor
{
        typedef FieldLines::IntegrationPoint<Eagle::PhysicalSpace::tvector> input_t;
        typedef Eagle::point3 output_t;

        output_t operator()(const input_t&in) const
        {
                return output_t(in.location);
        }
};

template <>
inline void CollectLineData(Representation&CartesianVertices, 
                            const std::vector<RefPtr<Chunk< FieldLines::IntegrationPoint<Eagle::PhysicalSpace::tvector> > > >&Lines, size_t nPoints)
{
RefPtr<Field>   LinePositions = new Field();
        LinePositions->setPersistentData( CollectLines( CopyExtractor(), Lines, nPoints) );
        CartesianVertices.setPositions( LinePositions );
}






struct  ExtractPointLocation
{
        typedef FieldLines::IntegrationPoint<Eagle::metric44> input_t;
        typedef Eagle::point3 output_t;

        output_t operator()(const input_t&in) const
        {
                return output_t(in.location[1],in.location[2],in.location[3]);
        }       
};

struct  ExtractTime
{
        typedef FieldLines::IntegrationPoint<Eagle::metric44> input_t;
        typedef double output_t;

        output_t operator()(const input_t&in) const
        {
                return in.location[0];
        }
};

template <>
inline void CollectLineData(Representation&CartesianVertices, 
                            const std::vector<RefPtr<Chunk< FieldLines::IntegrationPoint<Eagle::metric44> > > >&Lines, size_t nPoints)
{
RefPtr<Field>   LinePositions = new Field();
        LinePositions->setPersistentData( CollectLines( ExtractPointLocation(), Lines, nPoints) );
        CartesianVertices.setPositions( LinePositions );

        CartesianVertices[ "ProperTime" ]->setPersistentData( CollectLines( ExtractTime(), Lines, nPoints) );
}



template <class Container>
RefPtr<LineSet::LinesetArray_t> createLineSet(const std::vector<RefPtr<Chunk<Container> > >&Lines)
{
size_t  NumberOfLines = Lines.size();
Ref<LineSet::LinesetArray_t> LinesetArray(NumberOfLines);
MultiArray<1, std::vector<LineSet::LineIndex_t> >&LineIndices = *LinesetArray; 

index_t VertexNumber = 0;
        for(index_t line=0; line<NumberOfLines; line++)
        {
        std::vector<LineSet::LineIndex_t>&E = LineIndices[ line ];
        index_t LineLength = Lines[line]->std_vector().size();

                E.resize( LineLength );
                for(index_t i=0; i<LineLength; i++)
                {
                        E[i] = VertexNumber;
                        VertexNumber ++;
                }
        }

        return LinesetArray;
}




template <class FieldType, typename LineType, int InterpolType>
class   LineAdvancer
{
typedef RefPtr<Chunk<IntegrationPoint<FieldType> > > IntegrationResult_t;
typedef MemArray<3, FieldType> ToIntegrateArray_t; 

        AtomicIntegrator<FieldType, LineType, InterpolType> AI;
        std::vector<IntegrationResult_t> Lines;
        bool dop;
        size_t nPoints;
        double min_line_length;
        double start_time;
        double time;

public:
        LineAdvancer(const std::vector<IntegrationResult_t>&LinesP,
                     const FieldSelector&FieldSelectionP, const RefPtr<LocalFromWorldPoint>&LocalPointFinderP,
                     const RefPtr<BoundingBox>&BBP, const double step_sizeP, const bool dopP, const size_t nPointsP,
                     const double min_line_lengthP, const double start_timeP, const double scale_facP )
        : AI(LocalPointFinderP, FieldSelectionP, BBP, step_sizeP, scale_facP)
        , Lines(LinesP)
        , dop(dopP)
        , nPoints(nPointsP)
        , min_line_length(min_line_lengthP)
        , start_time(start_timeP)
        , time(start_timeP)
        {}

//      LineIntegrator::LineIntegrator(){}


        size_t getNPoints() { return nPoints;}
        
//      std::vector<IntegrationResult_t>& getLines(){ return Lines; }


        bool advanceBreadthFirst( )
        {
        unsigned end_counter = 0; 

                if( dop && !AI.inited )
                        AI.initDop853( Lines.size() ); 

                for(unsigned p = 0; p < Lines.size(); p++)
                {
//                      printf("\n WORKING ON LINE %d out of %d\n", p, Lines.size() ); 
//                      VActionNotifier::ProgressInfo("Working on line " + String(int(p)) + " out of " + String(int(Lines.size())) );

                IntegrationResult_t &IntegrationLine = Lines[ p ];
                int last_index = IntegrationLine->size()-1; 


                        if( last_index >= 0 )
                                if( (*IntegrationLine)[last_index].is_end == true )
                                        continue; 

                        if( last_index >= 1 )
                        {
                                //printf( "advacer: length %f, limit: %f\n", (*IntegrationLine)[last_index-1].length, min_line_length ); 
                                if( (*IntegrationLine)[last_index-1].length > min_line_length )
                                {
                                        (*IntegrationLine)[last_index].setEnd(); 
                                        continue;
                                }                       
                        }       

                bool test;

                        if( !dop )
                                test = AI.doEuler( *IntegrationLine, time); 
                        else
                                test = AI.doDop853( *IntegrationLine, time, p ); 

                        if( test == false ) 
                                end_counter++;
                        else
                                nPoints++; 

        } 

                VActionNotifier::ProgressInfo("Integration of " + String(int(Lines.size())) + " lines done."  );

//      printf("end_counter %d, Lines.size() %d\n", end_counter, Lines.size() ); 
                return ( end_counter == Lines.size() ) ? true : false;
        }

        void finalize()
        {
                for(unsigned p = 0; p < Lines.size(); p++)
                {
                IntegrationResult_t &IntegrationLine = Lines[ p ]; 

                        AI.finalize( *IntegrationLine, time ); 
                }
        }

};


// type trait to convert the integration line type to a string used for naming in the fiber bundle
// must be provided in the specializations
template<typename FieldType, typename LineType>
struct LineTypeString
{
        string name() { return "Unknown"; }
};


template <typename FieldType, typename LineType>
class   IntegralLines : public IntegralFace
{
public:
        typedef MemArray<3, double>                    ScalarArray_t;
        typedef MemArray<3, Eagle::point3>             CoordsArray_t;
        typedef MemArray<1, Eagle::point3>             CoordsArray1D_t;
        typedef MemArray<1, tvector>                   VecsArray1D_t;
        typedef MemArray<1, std::string>               FragmentIDs_t;
        typedef MemArray<1, BoundingBox>               FragmentBounds_t;


        IntegralLines(const string&name, int p, const RefPtr<VCreationPreferences>&VP)
        : VObject(name, p, VP)
        , Fish<VObject>(this)
        , Fish<Field>("inputfield")
        , IntegralFace(name, p, VP)
        {}

        typedef RefPtr<Chunk<IntegrationPoint<FieldType> > > IntegrationResult_t;

        typedef typename FieldType::Chart_t Chart_t;
        typedef typename Coordinates<Chart_t>::point point;     
        typedef typename Coordinates<Chart_t>::vector tvector4; 
        
        override bool update(VRequest&R, double precision);
        
        static  string createChildname(const string&inputfield)
        {
                LineTypeString<FieldType,LineType> tmp;
                return tmp.name() +"(" + inputfield + ")";
        }


        void printProgress( unsigned&old_nPoints, const unsigned nPoints, const unsigned nPointsEst, VRequest&Context)
        {
                if( old_nPoints + 100 < nPoints )
                {
                        old_nPoints = nPoints;
                 char buf[64]; 
                        float percent = (double)nPoints / nPointsEst * 100; 
                        sprintf(buf, "Integrating %.1f %% #%d", percent, (int)nPoints); 
                        VActionNotifier::ProgressInfo (string(buf)); 
                        printf("%s\n", buf);fflush(stdout);
                        setStatusInfo( Context, buf);
                 } 
                
        }

        template<class IntegratorType>
        int doMainLoop(IntegratorType&IntT, const double max_steps, const int nPointsEst, const size_t nPoints, VRequest&Context)
        {
        int step_counter = 0; 
        bool hasfinished = false; 
        unsigned old_nPoints = 0; 
        
                while( hasfinished == false && step_counter < max_steps )
                {
                        printProgress(old_nPoints, nPoints, nPointsEst, Context ); 
                        hasfinished = IntT.advanceBreadthFirst();//,                   //'returns' 
                        step_counter++; 
                } 
                IntT.finalize();//, 

                return IntT.getNPoints(); 
        }

};

template<typename FieldType, typename LineType>
bool IntegralLines<FieldType,LineType>::update(VRequest&Context, double precision)
{
double sec = 0.0;
VTimer Tim_update; 

#ifdef VERBOSE
        printf("IntegralLines::update()  -->%s<--: enter\n", Name().c_str() );fflush(stdout); 
#endif 

// 1. Get input field information 
// 
#ifdef VERBOSE
        puts("IntegralLines::update() 1. Get input field information");fflush(stdout); 
#endif

FieldSelector   FieldSelection = getFieldSelector(Context); 
RefPtr<Field>   InputToIntegrateField = FieldSelection.getField(); 

        if(!FieldSelection.theSourceBundle) { return errorMsg(Context, "No bundle found in integrate field");    }
        if(!InputToIntegrateField)          { return errorMsg(Context, "No integratable field in input bundle"); }

double time      = FieldSelection.getTime(); 
string fieldName = FieldSelection.getFieldName();


// 2. Get emitter info
// 

GridSelector EmitterGS; 
FieldSelector EmitterFS; 
        StartGrid << Context >> EmitterGS; 
        StartField << Context >> EmitterFS; 

EmitterFieldData EFD( EmitterGS, EmitterFS ); 

RefPtr<Field> EmitterCoordinates = getEmitterField( Context, EFD ); 

Enum Dop;
        UseDop853 << Context >> Dop; 

Enum Interpol; 
        Ipol << Context >> Interpol;
        

//
// 3. Construct name for the output grid (the integral lines) and check if it already exists. And if so, if not the 
//    input is newer, such that it would require re-creation.
// 
string grid_name = FieldSelection.Gridname() + "_" + Name(); 

        if (Dop("yes") ) grid_name += "_dop853";

GridSelector GS; 
        myIntegralLines << Context >> GS; 

UpdateCheckData UCD( EmitterCoordinates, FieldSelection, GS, InputToIntegrateField, time, grid_name); 

        if( !needUpdate( Context, UCD) )
        {
                puts("IntegralLines::update() 3.5 no update needed");
                return true;
        }
// 
// 4. Make sure we have some seed points 
// 
        puts("IntegralLines::update() 4. Make sure we have seeding points");fflush(stdout); 

RefPtr<BoundingBox> ToIntegratefieldBBox = getBoundingBox( *FieldSelection.getGrid() ); 
        if (!ToIntegratefieldBBox) return setStatusError(Context, "No bounding box available.\n"); 

int     NumberOfInternalSeedPoints = 23; 
        InternalEmissionSeedPoints << Context >> NumberOfInternalSeedPoints;
AutoEmitterData AED( EmitterCoordinates, ToIntegratefieldBBox, NumberOfInternalSeedPoints);
        setDefaultEmissionPoints(Context, AED);


        //
        // Note: the emitter coordinate field could also be fragmented.
        // In this case, the getData() call will return a NullPtr().
        // In general, we need to rather iterate here over all field
        // fragments of the emitter coordinates.
        //
RefPtr<CoordsArray1D_t > EmitterVerts = EmitterCoordinates->getData();
        if (!EmitterVerts)
                return setStatusError( Context, "Emitter points not in a supported format"); 

        MultiArray<1, Eagle::point3>&EmitterCrds = *EmitterVerts; 

//
RefPtr<Field> EmitterVecField = EmitterFS.getField(); 
RefPtr<VecsArray1D_t > EmitterVecs;
 
        if(EmitterVecField)
                EmitterVecs = EmitterVecField->getData(); 



// 

// 5. Prepare data structs for integration 
// 
// 
        puts("IntegralLines::update() 5. Prepare data");fflush(stdout); 
BundlePtr   outBPtr = FieldSelection.BundleSource();
Info<Slice> currentSlice = FieldSelection.FieldSource; 

        assert( FieldSelection.getGrid() );
const Grid&InputGrid = *FieldSelection.getGrid(); 

        setStatusInfo(Context, "Setting up index search"); 
double  UnimapperScale = 1.0; 
        UniMapperScalator << Context >> UnimapperScale;

RefPtr<LocalFromWorldPoint>
        LocalPointFinder = new LocalFromWorldPoint(*currentSlice.getSlice(), FieldSelection.getGrid(),
                                                   FieldSelection.getGridname(),
                                                   FieldSelection.getGrid()->getCartesianPositions(), 
                                                   UnimapperScale );

double  step_size;
double  min_line_length = 0.0; 
double  max_steps = 0.0; 
double  scale_fac = 1.0;
        MinLineLength  << Context >> min_line_length; 
        MaxSteps << Context >> max_steps;
        StepSize    << Context >> step_size; 
        ScaleFac    << Context >> scale_fac;

        step_size *= ToIntegratefieldBBox->radius();
        step_size /= 33.0;

Enum AddData, AddMagnitude; 
        IncludeData << Context >> AddData; 
        IncludeMagnitude << Context >> AddMagnitude; 

        if( max_steps < 2)
                max_steps = 2; 


unsigned int nEmPoints = EmitterCrds.Size()[0]; 

#ifdef VERBOSE
        cout << "IntegralLines::update() Number of Seeds: " <<  nEmPoints << endl; 
#endif

//index_t       NumberOfLines = nEmPoints;

// 
// 6. Do the actual integration
// 
        puts("IntegralLines::update() 6. Do the actual");fflush(stdout); 
        setStatusInfo(Context, "Integrating ...");

unsigned int nth =  nEmPoints/10; 
        nth = (nth == 0)?5:nth; 

size_t  nPoints = nEmPoints; // because start points are first pushed onto the line
 
//std::vector<IntegrationResult_t> Lines( nEmPoints );
 std::vector<IntegrationResult_t> Lines( nEmPoints );
//Lines.resize( nEmPoints ); 
 //printf("IntegralLines::update() Lines.size: %d\n", (int)Lines.size());fflush(stdout); 

// init Lines by pushing first point on line 
//for( unsigned i = 0; i < Lines.size(); i++ ) 
for( unsigned i = 0; i < Lines.size(); i++ )
{
        Lines[i] = new Chunk<IntegrationPoint<FieldType> >(0); 
        IntegrationPoint<FieldType> start;
 
        // awful handling of 4d and 3d points. 3d emitters copied to point4 T component set to 0
        for(unsigned ii = 0; ii < IntegrationPoint<FieldType>::point::Dims; ii++)
        {
                if( IntegrationPoint<FieldType>::point::Dims < 4 )
                        start.location[ii] = EmitterCrds[i][ii]; 
                else
                {
                        start.location[0]=0.0;
                        if( ii < 3)
                                start.location[ii+1] = EmitterCrds[i][ii];      
                } 
                
        } 
#ifdef VERBOSE  
        cout << start.location << endl; 
#endif
        if(EmitterVecs )
        {
                MultiArray<1, tvector>&EmitterCrds = *EmitterVecs; 

                for(unsigned ii = 0; ii < IntegrationPoint<FieldType>::point::Dims; ii++)
                {
                        if( IntegrationPoint<FieldType>::point::Dims < 4 )
                                start.direction[ii] = EmitterCrds[i][ii]; 
                        else
                        {
                                start.direction[0] = 0.0; // note that this is not correct for lightlike geodesiccs. 
                                                          // but since the metric has to be know, this is set to the 
                                                          // corect value in the geodesic44 integrator class
                                if( ii < 3)
                                        start.direction[ii+1] = EmitterCrds[i][ii];     
                        } 
                } 
        }
        start.is_end = false; 

        Lines[i]->push_back(start); 
//      cout << "IntegralLines::update() startpoint:"<< start.location << " with start dir:" << start.direction <<endl;
} 

puts("1"); 

typedef LineAdvancer<FieldType, LineType, FieldInterpolatorBase::LINEAR   > AdvancerLinear;  
typedef LineAdvancer<FieldType, LineType, FieldInterpolatorBase::CUBIC    > AdvancerCubic; 
typedef LineAdvancer<FieldType, LineType, FieldInterpolatorBase::ANALYTIC > AdvancerAnalytic; 

int nPointsEst = nEmPoints * ((int)max_steps); 

        // main loop
        puts("2");
        ________________________________________________getTime(sec , Tim_update ); 

        if( Interpol("linear") )
        {

        AdvancerLinear  MyAdvancerLinear( Lines, FieldSelection, LocalPointFinder, ToIntegratefieldBBox, 
                                          step_size, Dop("yes"), nPoints, min_line_length, time, scale_fac ); 

                nPoints = doMainLoop<AdvancerLinear>(MyAdvancerLinear, max_steps, nPointsEst, nPoints, Context); 
        }
        else if( Interpol("cubic") )
        {
        AdvancerCubic  MyAdvancerCubic( Lines, FieldSelection, LocalPointFinder, ToIntegratefieldBBox, 
                                        step_size, Dop("yes"), nPoints, min_line_length, time, scale_fac ); 

                nPoints = doMainLoop<AdvancerCubic>(MyAdvancerCubic, max_steps, nPointsEst, nPoints, Context); 
        } 
        else
        {
        AdvancerAnalytic  MyAdvancerAnalytic( Lines, FieldSelection, LocalPointFinder, ToIntegratefieldBBox, 
                                              step_size, Dop("yes"), nPoints, min_line_length, time, scale_fac ); 

                nPoints = doMainLoop<AdvancerAnalytic>(MyAdvancerAnalytic, max_steps, nPointsEst, nPoints, Context); 
        } 

        double overall_int_time = Tim_update.secs() - sec;       
        
//      ________________________________________________printTime("IntegralLines::update(): Compute Time:", Tim_update.secs() - sec); 

        printf("IntegralLines::update(): ComputeTime: %f\n",overall_int_time  ); 
        printf("IntegralLines::update(): Time per integration %d points: %f\n", (int)nPoints, overall_int_time/nPoints ); 

        CurviCellsPerUniCell << Context << LocalPointFinder->MaxListSize;


// 7. Collect all points into one array 
//    Alternatively, could store all data in field fragments in the 
//    output positions field. This were more efficient, as it does 
//    not require to copy the points here. However, the display methods 
//    such as line render, must be checked to be able to handle such 
//    line data. Nevertheless, this were the way how to do it "right". 
//    There should be another, generic, module anyway that combines 
//    fragmented fields into unfragmented ones, for those modules which
//    are incapable elsewhere to operate on fragmented fields.
// 
//      puts("IntegralLines::update() 7. Collect");fflush(stdout); 


// 
// 8. Bake a new Grid object that carries the actual lines 
// 
// 
        puts("IntegralLines::update() 8. Bake");fflush(stdout); 
        {
        //
        // a. Create a Grid object on the Bundle, compatible with the Bundle,
        //    but not yet inserted into the Bundle. It will thus be invisible
        //    until we are done with baking it. Otherwise, an incomplete Grid
        //    would be visible in the Bundle, and a concurrent access to the
        //    Bundle's Grid would fail (as it may happen in a multithreaded
        //    environment). Alternatively, one could start baking the Grid already
        //    during integration, and publish it in the Bundle. Then a concurrent
        //    rendering would show the active computation going on. However, the
        //    Grid must always be internally consistent. One must never publish
        //    a half-done Grid into the Bundle.
        //
RefPtr<Grid> IntegralLineGrid = FieldSelection.theSourceBundle->newGrid(); 


LineSet OutputLineSet( IntegralLineGrid,
                       NullPtr(), // coordinates will be added a bit later
                       createLineSet( Lines ) ); 

        assert(OutputLineSet.CartesianVertices);

        CollectLineData(*OutputLineSet.CartesianVertices, Lines, nPoints);

        //
        // d. Define relationship between the output grid of lines, and the 
        //    grid on which the vector field lives. This is optional, but will 
        //    easy further operations on the line grid, such as evaluating fields 
        //    from the vector field's grid along the lines.
        // 
                // Remember the source grid's fragment ID's for each point of the 
                // integration lines. This is useful for multiblock grids, such that 
                // subsequent lookup of a point's fragment ID is not necessary, for 
                // instance when another field of the source grid shall be evaluated 
                // on the lines. 
                if (0) // wieso geht des nit mit geodesics44?
                {
                RefPtr<Skeleton>SourceFragments = InputGrid( SkeletonID(3, 2) );
                        assert(SourceFragments);
                typedef MemArray<1, std::string> FragmentIDs_t;
                Ref<FragmentIDs_t> LinesVertFragID( nPoints );
                MultiArray<1,std::string > LineVFID = *LinesVertFragID;
                {
                index_t vert_count = 0;
                        for(index_t lin = 0; lin <  Lines.size(); lin++ )
                        {
                        const std::vector<IntegrationPoint<FieldType> >&Line = Lines[lin]->std_vector();
                                for(index_t ver = 0; ver < Line.size() ; ver++ )
                                {
                                        LineVFID[ MIndex(vert_count) ] =  Line[ver].fragment;
                                        vert_count++;
                                }
                        } 
                        //assert( vert_count == nPoints);
                } 

                Representation&LineVerticesAsSourceFragments = (*OutputLineSet.Vertices)[ SourceFragments ]; 
                        LineVerticesAsSourceFragments.setPositions( new Field( LinesVertFragID ) );
                }

        //
        // e. Copy floating point indices from the integration lines to the line grid.
        //    Same as with the fragment ID's, these indices correspond to a representation 
        //    of the output grid in terms of the input grid and will ease mapping fields 
        //    in a future operation.
        //
                {
                typedef MemArray<1, Eagle::point3> FloatIndices_t;
                Ref<FloatIndices_t> LinesVertFI( nPoints ); 
                MultiArray<1, Eagle::point3 > LineVFI = *LinesVertFI;
                {
                index_t vert_count = 0;
                        for(index_t lin = 0; lin < Lines.size(); lin++ )
                        {
                        const std::vector<IntegrationPoint<FieldType> >&Line = Lines[lin]->std_vector();

                                for(index_t ver = 0; ver < Line.size() ; ver++ )
                                {
                                        LineVFI[ MIndex(vert_count) ] =  Line[ver].index_location;
                                        vert_count++;
                                }
                        } 
                        //assert( vert_count == nPoints);
                } 

                Skeleton&SourceSkeleton                = *(FieldSelection.getGrid()->findVertices());
                Representation&LineVerticesAsSourceVertices = (*OutputLineSet.Vertices)[ SourceSkeleton ]; 
                        LineVerticesAsSourceVertices.setPositions(  new Field( LinesVertFI ) );

                // f. add direction (i.e. velocity)
                {
/*
  Look at the classes above on CollectLineData how to extract 3D information
  from nD point data IF required.


  The two if's down here seem to do exactly the same, except that in
  the 1st case it uses the point dims for the loop after verifying that
  Dims is three, in the 2nd case it writes 3 explicitly. What's the point
  of

  if (Dims==3) f(Dims);
  if (Dims==4) f(3);

  ??

 */

                        if(IntegrationPoint<FieldType>::point::Dims == 3)
                        {
                        Ref<MemArray<1, Eagle::tvector3> >  LinesData( nPoints ); 
                        MultiArray<1, Eagle::tvector3>&LineData = *LinesData; 
                        index_t vert_count = 0;
                                for(index_t lin = 0; lin <  Lines.size(); lin++ )
                                {
                                const std::vector<IntegrationPoint<FieldType> >&Line = Lines[lin]->std_vector(); 
                                        for(index_t ver = 0; ver < Line.size() ; ver++ )
                                        {
                                                for(unsigned j = 0; j < IntegrationPoint<FieldType>::point::Dims; j++)
                                                        LineData[ MIndex(vert_count) ][j] = Line[ver].direction[j]; 
                                                vert_count++;
                                        }
                                } 
                                //assert( vert_count == nPoints); 
                                (*OutputLineSet.CartesianVertices)[ LineSet::PredefinedFieldNames::Velocity ]->setPersistentData( LinesData );
                        } 
                        if(IntegrationPoint<FieldType>::point::Dims == 4)
                        {
                        // Ref<MemArray<1,Eagle::tvector4> >  LinesData( nPoints ); 
                        // MultiArray<1,Eagle::tvector4>&LineData = *LinesData; 
                        // index_t vert_count = 0; 
                        //      for(index_t lin = 0; lin <  Lines.size(); lin++ )
                        //      {
                        //      const std::vector<IntegrationPoint<FieldType> >&Line = Lines[lin]->std_vector(); 
                        //              for(index_t ver = 0; ver < Line.size() ; ver++ )
                        //              {
                        //                      for(unsigned j = 0; j < IntegrationPoint<FieldType>::point::Dims; j++)
                        //                              LineData[ MIndex(vert_count) ][j] = Line[ver].direction[j]; 
                        //                      vert_count++;
                        //              }
                        //      } 
                        //      assert( vert_count == nPoints); 
                        // string tmp = string("Direction");
                        //      VerticesAsCartesian[ tmp ]->setPersistentData( LinesData );
                        Ref<MemArray<1, Eagle::tvector3> >  LinesData( nPoints ); 
                        MultiArray<1, Eagle::tvector3>&LineData = *LinesData; 
                        index_t vert_count = 0; 
                                for(index_t lin = 0; lin <  Lines.size(); lin++ )
                                {
                                const std::vector<IntegrationPoint<FieldType> >&Line = Lines[lin]->std_vector(); 
                                        for(index_t ver = 0; ver < Line.size() ; ver++ )
                                        {
                                                for(unsigned j = 0; j < 3; j++)
                                                         LineData[ MIndex(vert_count) ][j] = Line[ver].direction[j+1]; 
                                                cout<< "InegralLine:update(): Direction: " << LineData[ MIndex(vert_count) ] << endl;   
                                                vert_count++; 
                                        }
                                } 
                                puts("directions done 1");fflush(stdout);
                                //assert( vert_count == nPoints); 
                                (*OutputLineSet.CartesianVertices)[  LineSet::PredefinedFieldNames::Velocity ]->setPersistentData( LinesData ); 
                                puts("directions done 2");fflush(stdout);
                        }       
                        }
                } 

                
        // 
        // g. Add the data field to the output grid of lines. // this maybe needs special treatment 
        //
                if( AddData("yes") )
                {
                typedef MemArray<1, FieldType> LinesDataArr_t; 
                Ref<LinesDataArr_t>  LinesData( nPoints ); 
                MultiArray<1,FieldType>&LineData = *LinesData; 
                index_t vert_count = 0; 
                        for(index_t lin = 0; lin <  Lines.size(); lin++ )
                        {
                        const std::vector<IntegrationPoint<FieldType> >&Line = Lines[lin]->std_vector(); 
                                for(index_t ver = 0; ver < Line.size() ; ver++ )
                                {
                                        LineData[ MIndex(vert_count) ] = Line[ver].data;
                                        vert_count++;
                                }
                        } 
                        //assert( vert_count == nPoints);
                        (*OutputLineSet.CartesianVertices)[ Fieldname(Context) ]->setPersistentData( LinesData );
                } 



                        // And establish the same vector field as the tangential vectors 
                        // given on this line set. This is an alias to the same vector field, 
                        // under which TangentialVectors operations will find this field by 
                        // convention as defined in fiber/baseop/TangentialVectors.hpp 
                        // 
                        //VerticesAsCartesian[ TangentialVectorFieldName ]->setPersistentData( LinesData ); 

        // 
        // h. As a feature of this module, compute the magnitude of the vector field // needs special treatment
        //    and store it in the lines. 
        // 
                if (AddMagnitude("yes")) 
                {
                Ref<MemArray<1, double> > LinesVectorMagnitude(nPoints); 
                MultiArray<1,double >&LineMag = *LinesVectorMagnitude; 
                index_t vert_count = 0;
                        for(index_t lin = 0; lin < Lines.size(); lin++ )
                        {
                        const std::vector<IntegrationPoint<FieldType> >&Line = Lines[lin]->std_vector(); 
                                for(index_t ver = 0; ver < Line.size() ; ver++ )
                                {
                                        LineMag[ MIndex( vert_count) ] = norm( Line[ver].data );
                                        vert_count++; 
                                }
                        }

#define DOUBLELINE "\u2016"
                string outFieldName( DOUBLELINE + Fieldname(Context) + DOUBLELINE ); 
                        (*OutputLineSet.CartesianVertices)[ outFieldName ]->setPersistentData( LinesVectorMagnitude );
//              puts("IntegralLine::update() DEB: setPersistent LineVectorMagnitudesF");
                } 

        // 
        // y. Allow standard fields from the LineSet class to be available on integrated data
        // 
                puts("before setup standard");fflush(stdout);
                if (!OutputLineSet.setupStandardFields( IntegralLineGrid ) )
                {
                        puts("LINE INTEGRATOR ERROR: Could not set up standard fields on lines!");fflush(stdout);
                }



        // 
        // z. Finally, insert the Grid object into the bundle, at the current 
        //    time slice. It will hereby be visible in the bundle and ready 
        //    to be used elsewhere (like a rendering module). At last, announce 
        //    this Grid's properties as the output of this VISH Object.
        // 

                puts("insert slice");fflush(stdout);
                currentSlice.getSlice()->insert( grid_name, IntegralLineGrid ); 

                FieldSelection.BundleSource().speak("FieldSelection.BundleSource()"); 

                
        GridSelector GS( grid_name, FieldSelection.BundleSource() ); 
                myIntegralLines << Context << GS;
        } 

        infoMsg(Context, "Done, output at: " + grid_name);

//      printf("ComputMultipleStreamLines::update() Written into Bundle, setPersitentData and touched Spacetimeslot()\n");
        VActionNotifier::ProgressInfo (string("Streamline computation complete")); 

        puts("IntegralLine::update() >>>  EXIT");

        return true;
}

}