IntegralHeart.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 <fish/fiber/vector/Interpolate.hpp>
#include <fish/fiber/vector/LinearIpol.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 <fish/fiber/baseop/CopySkeletons.hpp>
//#include <fish/fiber/baseop/AMRFieldFinder.hpp>

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

#include <ocean/plankton/VTimer.hpp>

#include <grid/types/LineSet.hpp>

//#define VERBOSE

#define TIMEIT


namespace Fiber
{
        // convinience function, to be placed somewhere else
        template<int N, typename T>
        RefPtr<MemArray<N, T> > getMemArrayViaMBase( RefPtr<Field>& myField, RefPtr<FragmentID> FragID = MemCore::NullPtr() )
        {
                RefPtr<CreativeArrayBase> cBase = myField->getCreator(); 
                if(!cBase)
                {
                        puts("getMemArrayViaMBase() No cBase for FIELD found in field");fflush(stdout); 
                        return MemCore::NullPtr(); 
                }       

                RefPtr<MemBase> mBase = cBase->create(); 
                if(!mBase)
                {
                        puts("getMemArrayViaMBase() No mBase for FIELD found in field");fflush(stdout); 
                        return MemCore::NullPtr(); 
                } 

                RefPtr<MemArray<N, T> > ret = mBase; 
                ret = myField->getCreator(FragID)->create(); 
                if(!ret)
                {
                                mBase.speak("MBase Inforation ------>"); 
                                puts("getMemArrayViaMBase() ERROR retrieving matching field."); 

                                return MemCore::NullPtr();
                } 

                return ret;
        }

        template<typename T>
        RefPtr<MemCore::TypedChunk<T> > getChunkOfField( RefPtr<Field>& myField, RefPtr<FragmentID> FragID = MemCore::NullPtr() )
        {
                RefPtr<MemArray<1, T> > DataArr = getMemArrayViaMBase<1, T>( myField ); 
                assert(DataArr); 

                return DataArr->myChunk();
        }
}

namespace FieldLines
{
typedef std::string FragmentID_t;

// This class iterates over the input grid object and
// maps Field with time value for time dependent fields.
//template<typename FieldType, typename LineType>
class GridIterator : public EvolutionIterator<Grid>
                   , public MemCore::ReferenceBase<GridIterator>
{
        typedef std::map<double, RefPtr<Field> > MapFieldObject_t;      
        std::string theFieldName;
public:
        MapFieldObject_t fieldObjectCollection;
        std::map<double, double> timeDifferenceMap;
        double t0, t1;

        // Constructor
        GridIterator(const string fieldName)
        : MemCore::ReferenceBase<GridIterator>(this)
        , theFieldName(fieldName)
        , t0(0.0)
        , t1(0.0)
        {}
        
        // Define the apply function here.
        bool apply(double t, GridID&id, Grid&g)
        {       
                t1 = t;
                if(t1 != t0)
                {
                        double stepTime = t1-t0;
                        timeDifferenceMap[t1] = timeDifferenceMap[t0] = stepTime;
                }
                else 
                        timeDifferenceMap[t0] = 0.0;
                        
                t0 = t;
                
                // Get the Field object with corresponding fieldName
                RefPtr<Field> fieldObj = g(theFieldName); 

                if(!fieldObj)
                {
                        printf("GridIterator::apply() Null to field retrieved. Abording at t=%f!", t); 
                        return true;
                }

//              assert(fieldObj);
                
                fieldObjectCollection[t] = fieldObj; 

#ifdef VERBOSE
                printf("GridIterator::apply() inserted %s st %.18f --> map size:%d\n", theFieldName.c_str(), t, (int)fieldObjectCollection.size() ); 
                fflush(stdout);
#endif
                return true;    
        }               
};


class FieldCollectionDataBase1D
{
protected:
        RefPtr<Grid> myGrid;
        RefPtr<MemArray<1, Eagle::point3> > Positions;

        unsigned size;

//      RefPtr<MemArray<1, Eagle::point3> > LocalPositions;
//      RefPtr<MemArray<1, std::string> >   FragmentNames;
public:
        FieldCollectionDataBase1D( const RefPtr<Grid>& myGridP )
        : myGrid(myGridP)
        , size(0)
        {
                RefPtr<Field> PosField = myGrid->getCartesianPositions(); 

                if( !PosField ) 
                        puts("FieldCollectionDataBase1D::FieldCollectionDataBase1D() ERROR No Field in Grid"); 
                else
                {
                        Positions = getMemArrayViaMBase<1, Eagle::point3>( PosField ); 

                        if( Positions )
                                size = Positions->Size()[0];
                }
        }

};

template <class FieldType, typename LineType>
struct FieldCollection : public MemCore::ReferenceBase<FieldCollection<FieldType,LineType> >
{
        FieldCollection( const RefPtr<Grid>& myFieldP)
        : MemCore::ReferenceBase<FieldCollection<FieldType,LineType> >(this)
        {}
        
        unsigned getVertexSize()
        {
                puts("FieldCollection::getVertexSize() not implemented, returning ZERO"); 
                return 0;
        }
        
        unsigned getDataSize()
        {
                puts("FieldCollection::getDataSize() not implemented, returning ZERO"); 
                return 0;
        }
};

struct AtomicDataBase
{
        RefPtr<LocalFromWorldPoint> LocalPointFinder;
        FieldSelector FieldSelection;
        RefPtr<BoundingBox> BB;
        RefPtr<GridIterator> GI;
        double step_size; 
        double scale;
        double nextTime;

        AtomicDataBase( const RefPtr<LocalFromWorldPoint>&LocalPointFinderP, 
                        const FieldSelector&FieldSelectionP, const RefPtr<BoundingBox>&BBP, const RefPtr<GridIterator>&GridIteratorP,
                        const double step_sizeP , const double scale_facP )
        : LocalPointFinder(LocalPointFinderP)
        , FieldSelection(FieldSelectionP)
        , BB(BBP)
        , GI(GridIteratorP)
        , step_size(step_sizeP)
        , scale(scale_facP)
        {}

        ~AtomicDataBase()
        {
#ifdef VERBOSE          
                puts("~AtomicDataBase()");fflush(stdout);
                LocalPointFinder.speak("--------------------------------------------------------------~AtomicDataBase() LocalPointFinder"); fflush(stdout);
                BB.speak("--------------------------------------------------------------~AtomicDataBase() BB"); fflush(stdout);
                GI.speak("--------------------------------------------------------------~AtomicDataBase() GI"); fflush(stdout);
#endif
        }
};

// type trait for atomic integration requires Euler and Dop853 functions
// the interpolation type is defined at compile time via the last template paramter.
template<typename FieldType, typename LineType, int InterpolType>
class AtomicIntegrator
{
        typedef MemArray<3, FieldType> ToIntegrateArray_t; 
public:
        
        virtual ~AtomicIntegrator()
        {
#ifdef VERBOSE  
                puts("~AtomicIntegrator() base");fflush(stdout);
#endif
        }

        unsigned long nextStartIndexPerCoarseCall( unsigned long start, unsigned long size )
        {
                puts("FieldLines::AtomicIntegrator<FieldType>.nextStartIndexPerCoarseCall() ERROR: Unsupported type Field/Line Type doing nothing!");fflush(stdout); 
                return 0;
        }

        // does the local data extraction from using the positions field of the front to comput local coords, fragmentIds
        // and interpolated field data.
        bool extractLocalData( FieldCollection<FieldType, LineType>& CurrentFront,
                               const unsigned index, const double time  )
        {
                puts("FieldLines::AtomicIntegrator<FieldType>.extractLocalData() ERROR: Unsupported type Field/Line Type doing nothing!");fflush(stdout); 
                return true;
        }

        bool doEuler( FieldCollection<FieldType, LineType>& lastFront,
                      unsigned index, 
                      FieldCollection<FieldType, LineType>& newFront,
                      double time,
                      unsigned integration_width)
        {
                puts("FieldLines::AtomicIntegrator<FieldType>.doEuler (Fronts) () ERROR: Unsupported type Field/Line Type doing nothing!");fflush(stdout); 
                return false;   
        }

        bool doDop853( FieldCollection<FieldType, LineType>& lastFront,
                       unsigned index, 
                       FieldCollection<FieldType, LineType>& newFront,
                       double time,
                       unsigned integration_width )
        {
                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;
        }       
        
        double getNextTime()
        {
                return 0.0;
        }
};

        
template <class FieldType, typename LineType>
class GridOperator
{
        RefPtr<Bundle>& myBundle;
        
public:
        GridOperator( RefPtr<Bundle> myBundleP = MemCore::NullPtr() )
        :myBundle( myBundleP )
        {}

        RefPtr<Grid> getPreparedGrid( const RefPtr<Grid>& EmitterGrid, unsigned long steps_estimate = 500 )
        {
                puts("GridOperator::getpreparedGrid() WARNING: Using default base implementation, returning new Grid"); 
                return myBundle->newGrid();
        }

        RefPtr<Grid> advanceGrid(const RefPtr<Grid>& gridP )
        {
                puts("GridOperator::advanceGrid() WARNING: Using default implementation, returning new Grid"); 
                return myBundle->newGrid();
        }

        RefPtr<Grid> refineGrid(const RefPtr<Grid>& gridP )
        {
                puts("GridOperator::advanceGrid() WARNING: Using default implementation, returning new Grid"); 
                return gridP;
        }

        void storeGrid( const std::string& grid_name, const RefPtr<Grid>& gridP, double time )
        {
                puts("GridOperator::storeGrid() WARNING: Using default implementation, doing Nothing"); 
        }

        void finalizeGrid( const std::string& grid_name, const RefPtr<Grid>& gridP, double slice )
        {
                puts("GridOperator::finalizeGrid() WARNING: Using default implementation, doing Nothing"); 
        }
        
        unsigned getIntegrationSize(  const RefPtr<Grid>& gridP )
        {
                puts("GridOperator::getIntegrationSize() WARNING: Using default implementation, returning ZERO"); 
                return 0;
        }

        RefPtr<Bundle> getBundle()
        {
                return myBundle;
        }
};



template <class FieldType, typename LineType, int InterpolType>
class  CoarseIntegrator
{
typedef MemArray<3, FieldType> ToIntegrateArray_t; 

        AtomicIntegrator<FieldType, LineType, InterpolType> AI;

        bool dop;
        size_t nPoints;
        double line_length;
        double start_time;
        double time;
        int number_lines;

public:
        CoarseIntegrator( const FieldSelector&FieldSelectionP, const RefPtr<LocalFromWorldPoint>&LocalPointFinderP,
                       const RefPtr<BoundingBox>&BBP, const RefPtr<GridIterator>&GridIteratorP, 
                       const double step_sizeP, const bool dopP, const size_t nPointsP,
                       const double line_lengthP, const double start_timeP, const double scale_facP )
        : AI( LocalPointFinderP, FieldSelectionP, BBP, GridIteratorP, step_sizeP, scale_facP )
        , dop(dopP)
        , nPoints(nPointsP)
        , line_length(line_lengthP)
        , start_time(start_timeP)
        , time(start_timeP)
        {}

        virtual ~CoarseIntegrator()
        {
#ifdef VERBOSE  
                puts("~FronAdvancer()");fflush(stdout);
#endif
        }

//      LineIntegrator::LineIntegrator(){}


        size_t getNPoints() { return nPoints;}

        bool advance( RefPtr<Grid>&CurrentGrid, RefPtr<Grid>&NewGrid,
                      double&time, unsigned long&start, unsigned long&size)
        {
//              puts("CoarseIntegrator::advance() enter" );fflush(stdout);
                FieldCollection<FieldType, LineType> CurrentFields( CurrentGrid ); 
//              puts("CoarseIntegrator::advance() 1" );fflush(stdout);
                FieldCollection<FieldType, LineType> NewFields( NewGrid ); 
//              puts("CoarseIntegrator::advance() 2" );fflush(stdout); 
                
                unsigned end_counter = 0; 
                
                if(dop)
                        AI.initDop853( size ); 
                        
                //printf("CoarseIntegrator::Time at the start of advnace = %.18f\n", time); 
                //printf("About to get into Euler...start = %d, size = %d\n", int(start), int(size) );

                // add a omp pragma here maybe
                for(unsigned p = start; p < start + size; p++)
                {
#ifdef VERBOSE
                        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(size)) ); 
                        printf("Working on point %d, about to get into Euler...start = %d, size = %d\n", p, int(start), int(size)); 
#endif
                        
                        bool test; 
                                
                        // could this condition also be pulled out to one condition in the update, like the integration type?
                        if( !dop )
                                test = AI.doEuler( CurrentFields, p, NewFields, time, size ); 
                        else
                                //      test = AI.doDop853( CurrentFields, p, NewFields, time, size ); 


                        if( test == false ) 
                                end_counter++;
                        else
                                nPoints++;
                } 
                // change the time value here ...MAY BE!!!
                time = AI.getNextTime();
                
#ifdef VERBOSE  
                printf("CoarseIntegrator::Time at the end of advance = %.18f\n", time);
#endif          
                // change the start index in the needs of the Atomic Integrator
                start = AI.nextStartIndexPerCoarseCall( start, size );

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

//      printf("end_counter %d, Lines.size() %d\n", end_counter, Lines.size() ); 

//      bool check =  ( end_counter == size ) ? true : false; 

//      printf("CoarseIntegrator::advance() exit, check: %d\n", check);fflush(stdout); 

                return true;
        }

        // should be implementable by filling the data arrays up to the size of the vertices
        bool extractLocalGridData( RefPtr<Grid > & CurrentGrid, double time )
        {
#ifdef VERBOSE  
                puts("CoarseIntegrator::extractLocalGridData() enter");fflush(stdout);
#endif
                FieldCollection<FieldType, LineType>CurrentFields( CurrentGrid ); 

#ifdef VERBOSE  
                puts("CoarseIntegrator::extractLocalGridData() 1");fflush(stdout);
#endif

                bool check = true; 

                for( unsigned p = CurrentFields.getDataSize(); p < CurrentFields.getVertexSize(); p++ )
                {
                        check = check && AI.extractLocalData( CurrentFields, p, time );
                }
 
#ifdef VERBOSE  
                printf("CoarseIntegrator::extractLocalGridData() exit, check: %d\n", check);fflush(stdout);
#endif

                return check;
        }

};


// 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   IntegralHeart : 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, Eagle::tvector3>           VecsArray1D_t;
        typedef MemArray<1, std::string>               FragmentIDs_t;
        typedef MemArray<1, BoundingBox>               FragmentBounds_t;
        
        struct FieldState : State
        {
                WeakPtr<GridIterator> flowIterator;
                double time;
                
                FieldState()
                : time(0.0)
                {}
        };
        
        override RefPtr<State> newState() const
        {
                return new FieldState();
        }


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

        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);

                
        struct  AMRLevelExtractor : public LevelIterator
        {
                double          currenttime;
                std::string     AMRfieldname;
                
                // store available refinement levels of Field in map (at current time)
//              std::map<unsigned, MemCore::RefPtr<Fiber::Field> > collect;
                std::map<double, std::map<unsigned, RefPtr<Fiber::Field> > > collect;

                AMRLevelExtractor(std::string& fieldnameP, double currenttimeP, std::map<double, std::map<unsigned, RefPtr<Fiber::Field> > >& collectP)
                : currenttime(currenttimeP)
                , AMRfieldname(fieldnameP)
                , collect(collectP)  
                {}

                // collect time, levelnr and field in map
                override bool apply(const RefPtr<Fiber::Representation>&CartesianLevelRep,
                                    double time, int Level,
                                    const RefPtr<Fiber::Grid>&GridWithCoarsetRefinementLevel, 
                                    const RefPtr<ValuePool>&Context)
                {
                RefPtr<Field> Coords     = CartesianLevelRep->Positions(); 
                RefPtr<Field> DataField  = (*CartesianLevelRep)(AMRfieldname); 

                        collect.insert( pair<unsigned, RefPtr<Fiber::Field> >( Level, DataField ) );
                        
                        return true;
                }

        }; 
        

        class AMRFieldFinder
        {
                Fiber::BundlePtr bundle;
                std::map<double, std::map<unsigned, RefPtr<Fiber::Field> > > AMRField;
                std::string fieldname;
                
                // iterate over grids in time
                struct AMRGridEvolution : EvolutionIterator<Fiber::Grid>
                {
                        std::string fieldname;
//                      std::map<double, std::map<unsigned, RefPtr<Fiber::Field> > > myAMRField;

                        AMRGridEvolution( std::string& fieldnameP) //, std::map<double, std::map<unsigned, RefPtr<Fiber::Field> > >&myAMRFieldP )
                        : fieldname(fieldnameP)
//                      , myAMRField(myAMRFieldP)
                                {};
                        
                        // first collect level data in map by using the level iterator. then store available field ptr and level nr 
                        // in time-map of Field finder
                        override bool apply (double t, GridID&id, Grid&g)
                                {
                                        std::map<unsigned, MemCore::RefPtr<Fiber::Field> > collect; 

                                        AMRLevelExtractor refIt(fieldname, t, collect); 

                                                g.iterate( refIt ); 

                                                AMRField.insert( pair<double, RefPtr<Field> >(t, collect) );
                                        
                                        return true;
                                }
                };
                
                
                
        public:
        
                AMRFieldFinder(const string& fieldnameP)
                : fieldname(fieldnameP)
                        {}
                
                // collect all time and field data into a the map structure
                void precompute()
                        {
                        }
                
                /* also provide a function to find field on demand without pre-iterations
                   pre evaluation might be not too efficient in huge time dependent data sets
                RefPtr<Fiber::Field> getFinestField(double time)
                {
                return MemCore::NullPtr();
                }
                */
                
        };
        
        static  string createChildname(const string&inputfield)
        {
                LineTypeString<FieldType,LineType> tmp;
                return tmp.name() +"(" + inputfield + ")";
        }

        // print progress of iteration as status and stdout
        void printProgress( unsigned long&old_nPoints, const unsigned long nPoints, const unsigned long 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("IntegralHeart::printProgress(): %s\n", buf);fflush(stdout);
                        setStatusInfo( Context, buf); 
                 } 
                
        }

        
        // Function to simplify main iteration loop in update.
        template<class IntegratorType>
        int doMainLoop(IntegratorType&IntT, 
                       RefPtr<Grid>& StartGrid,
                       GridOperator<FieldType, LineType>& GridOperator,
                       const string& grid_name, 
                       FieldSelector& FieldSelection,
                       double time, 
                       unsigned long max_steps, int nPointsEst, size_t nPoints, VRequest&Context)
        {
//              puts("IntegralHeart::doMainLoop() enter");fflush(stdout);
                
        unsigned long step_counter = 0; 
        unsigned long integration_width = GridOperator.getIntegrationSize( StartGrid ); 
        unsigned long start_count = 0;
        bool hasfinished = false; 
        unsigned long old_nPoints = 0; 

        RefPtr<Grid> CurrentGrid = StartGrid; 
        RefPtr<Grid> NewGrid;

                while( hasfinished == false && step_counter < max_steps )
                {
                        printProgress(old_nPoints, nPoints, nPointsEst, Context ); 
#ifdef VERBOSE  
                        puts("fill local data and store grid  in case of being the initial front");
#endif
                        if( step_counter == 0)
                        {                 
                                printf("Emitter Grid:: Time at which emitter grid is stored %.18f\n", time);                    
                                hasfinished = hasfinished || 
                                              !IntT.extractLocalGridData( CurrentGrid, time ); 
                                GridOperator.storeGrid( grid_name, CurrentGrid, time );
                        } 
#ifdef VERBOSE  
                        puts("maybe create a copy for the new grid");
#endif
                        NewGrid = GridOperator.advanceGrid( CurrentGrid ); 
#ifdef VERBOSE  
                        puts("do one integration step of the whole front"); 
#endif
                        hasfinished = hasfinished || 
                                      !IntT.advance( CurrentGrid, NewGrid, time, start_count, integration_width ); 
#ifdef VERBOSE                          
                        puts("refine the grid");
#endif
                        NewGrid = GridOperator.refineGrid( NewGrid ); 

#ifdef VERBOSE
                        puts("fill local data in front"); 
#endif                  
                        hasfinished = hasfinished || 
                                      !IntT.extractLocalGridData( NewGrid, time ); 
#ifdef VERBOSE                  
                        puts("store grid");
#endif
                        GridOperator.storeGrid( grid_name, NewGrid, time ); 

                        CurrentGrid = NewGrid; 

                        step_counter++; 
#ifdef VERBOSE
                        printf("hasfinished: %d, step_counter: %ld, max_steps: %ld \n\n", hasfinished, step_counter,  max_steps); fflush(stdout);
#endif
                } 

#ifdef VERBOSE
                puts("do a final operation on the grid");
#endif
                GridOperator.finalizeGrid( grid_name, NewGrid, time ); 

#ifdef VERBOSE
                puts("exit main loop");
#endif
                return IntT.getNPoints(); 
        }

};

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

        printf("IntegralLines::update()  -->%s<--: enter\n", Name().c_str() );fflush(stdout); 
 
// 1. Get input field information 
// 
        puts("IntegralLines::update() 1. Get input field information");fflush(stdout); 
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();

// Get the State Variable
RefPtr<FieldState> state = newState();


// 2. Get emitter info
// 
GridSelector EmitterGS; 
FieldSelector EmitterFS; 
        StartGrid << Context >> EmitterGS; 
        StartField << Context >> EmitterFS; 

        EmitterFieldData EFD( EmitterGS, EmitterFS ); 

        RefPtr<Field> EmitterCoordinates = getEmitterField( Context, EFD );
        
//
// 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(); 

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   bundlePtr = FieldSelection.BundleSource(); // Was named as outBPtr. Modified by Bidur
Info<Slice> currentSlice = FieldSelection.FieldSource; 

        assert( 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; 

        if( max_steps < 2)         max_steps = 2; 
        if( min_line_length < 0.0) min_line_length = 0.0; 

unsigned long max_steps_l = (unsigned long)std::floor(max_steps);

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

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

        printf("IntegralHeart::update() nEmPoints: %d, nPointsEst: %d\n",nEmPoints,  nPointsEst);

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

size_t nPoints = 0; 

        GridOperator<FieldType, LineType> integrationGridOperator( FieldSelection.getBundle() ); 

// TODO make sure emitter grid is there. overrides current version of emitter-grid-point-creation! 
RefPtr<Grid> IntegrationGeometry = integrationGridOperator.getPreparedGrid( EmitterGS(time, EmitterGS.getBundle() ), (unsigned long) max_steps );

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

// Define new GridIterator and assign it the Iterator from State variable 
RefPtr<GridIterator> flowIterator = state->flowIterator;

int totalSteps;
double gridTime; 

        // Declare and initialize Grid Iterator 
        flowIterator = new GridIterator( fieldName );
        gridTime = time;

        // This iterates over Grid objects of input Bundle.
        // FlowIterator then maps the time of Grid with the Field(vector field) object of Grid. 
        totalSteps = bundlePtr->iterateForward(gridTime, FieldSelection.Gridname(), *flowIterator, false); 
        printf("\nTotal time slices = %d and Time value = %f %f\n", totalSteps, time, gridTime); 


Enum Interpol; 
        Ipol << Context >> Interpol; 
Enum Dop;
        UseDop853 << Context >> Dop; 

        puts("2");

        ________________________________________________getTime(sec , Tim_update ); 

        if( Interpol("linear") ) 
        {

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

                nPoints = doMainLoop<AdvancerLinear>( MyAdvancerLinear, IntegrationGeometry, integrationGridOperator, grid_name,
                                                      FieldSelection, time, max_steps_l, nPointsEst, nPoints, Context); 
        
                nPoints = MyAdvancerLinear.getNPoints(); 
        } 
        else if( Interpol("cubic") )
        {
        AdvancerCubic  MyAdvancerCubic( FieldSelection, LocalPointFinder, ToIntegratefieldBBox, flowIterator,
                                          step_size, Dop("yes"), nPoints, min_line_length, time, scale_fac ); 

                nPoints = doMainLoop<AdvancerCubic>( MyAdvancerCubic, IntegrationGeometry, integrationGridOperator, grid_name, 
                                                     FieldSelection, time, max_steps_l, nPointsEst, nPoints, Context); 
                nPoints = MyAdvancerCubic.getNPoints(); 
        } 
        else
        {
        AdvancerAnalytic  MyAdvancerAnalytic( FieldSelection, LocalPointFinder, ToIntegratefieldBBox, flowIterator,
                                              step_size, Dop("yes"), nPoints, min_line_length, time, scale_fac ); 

                nPoints = doMainLoop<AdvancerAnalytic>( MyAdvancerAnalytic, IntegrationGeometry, integrationGridOperator, grid_name, 
                                                        FieldSelection, time, max_steps_l, nPointsEst, nPoints, Context ); 
                nPoints = MyAdvancerAnalytic.getNPoints(); 
        } 

#ifdef VERBOSE  
        puts("3");fflush(stdout);
#endif

        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)nPointsEst, overall_int_time / nPointsEst ); 

        CurviCellsPerUniCell << Context << LocalPointFinder->MaxListSize; 

#ifdef VERBOSE  
        puts("4");fflush(stdout);
#endif

        {
        GridSelector GS( grid_name, integrationGridOperator.getBundle() ); 
                myIntegralLines << Context << GS; 

                printf("time: %f\n", time);

/*
        Grid& myGrid = *GS(time); 

        puts("5");fflush(stdout);

        RefPtr<Field> myField = myGrid.CartesianPositions(); 

        if(!myField)
        {
                puts("No FIELD!"); return true;
        }

        RefPtr<MemArray<1, Eagle::point3> > myArr = getMemArrayViaMBase<1, Eagle::point3>( myField ); 
        if(!myArr)
        {
                puts("No ARR!"); return true;
        }

        MultiArray<1, Eagle::point3>&multArr = *myArr; 

        

        MultiIndex<1> mx = multArr.Size(); 
        printf("MultiArray: Size(): %ld\n", mx[0] ); 
        
        MultiIndex<1> mi; 

        unsigned long hyperlen = multArr.getLength(); 
        printf("MultiArray: hyper size: %ld\n", hyperlen ); 

        unsigned long hypercount = multArr.maxcount(); 
        printf("MultiArray: hyper size: %ld\n", hyperlen ); 
        printf("MultiArray: hyper count: %ld\n", hypercount ); 

        printf("MultiArray: multidim size: %ld\n", mx[0]); 

        for(unsigned i = 0; i < mx[0]; i++)
        {
                mi[0] = i;
                std::cout << multArr[mi] << endl;
        }
        printf("ArraySize: %ld\n", mx[0]); 

        RefPtr<Chunk<Eagle::point3> >PositionsChunk = getChunkOfField<Eagle::point3>( myField ); 
        if(!PositionsChunk)
        {
                puts("No CHUNK!"); return true;
        }

                for(unsigned i = 0; i < PositionsChunk->size(); i++)
                        std::cout << (*PositionsChunk)[i] << endl;
*/
        }

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

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


        return true;
}

}