IntegralHeart.hpp

#define ENABLE_MEMCORE_GETTIME
#include <memcore/Timer.hpp>
 
#include "IntegralFace.hpp"
 
#include <BaseSpace/FieldInterpolator.hpp>
#include <eagle/PhysicalSpace.hpp>
#include <eagle/STA.hpp>
#include <ocean/shrimp/VEnum.hpp>
#include <ocean/shrimp/PhysicalSpace.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 <memcore/Timer.hpp>
 
#include <grid/types/LineSet.hpp>
#include <fish/fiber/baseop/RefineSurface.hpp>
 
#ifdef  USE_OPENMP
#include <omp.h>
#endif
 
//#define VERBOSE
#define PRT_TIME_ONLY
#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 nextTime;
        unsigned long stepNr;
 
        AtomicDataBase( const RefPtr<LocalFromWorldPoint>&LocalPointFinderP,
                        const FieldSelector&FieldSelectionP, const RefPtr<BoundingBox>&BBP, const RefPtr<GridIterator>&GridIteratorP,
                        const double step_sizeP)
        : LocalPointFinder(LocalPointFinderP)
        , FieldSelection(FieldSelectionP)
        , BB(BBP)
        , GI(GridIteratorP)
        , step_size(step_sizeP)
        , nextTime(0.0)
        , stepNr(0)
        {}
 
        ~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>
struct GridOperatorData : public MemCore::ReferenceBase<GridOperatorData<FieldType,LineType> >
{
        GridOperatorData()
        : MemCore::ReferenceBase<GridOperatorData<FieldType,LineType> >(this)
         {}
 
};
 
template <class FieldType, typename LineType>
class GridOperator
{
        RefPtr<Bundle>& myBundle;
        RefPtr< GridOperatorData<FieldType, LineType> >& myData;
 
public:
//      GridOperator( RefPtr<Bundle> myBundleP = MemCore::NullPtr() )
//      :myBundle( myBundleP )
//      {}
 
        GridOperator( RefPtr<Bundle> myBundleP, RefPtr< GridOperatorData<FieldType, LineType> >& myDataP )
        :myBundle( myBundleP )
        ,myData( myDataP )
        {}
 
        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 )
        : AI( LocalPointFinderP, FieldSelectionP, BBP, GridIteratorP, step_sizeP )
        , 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
// #ifdef USE_OPENMP
 
// #ifdef VERBOSE
 
//              const int max_threads_possible =  omp_get_max_threads();
//              cout << "advance() using threads per front: " <<  max_threads_possible << endl;
 
//              exit(0);
 
// #endif
 
//      #pragma omp parallel
//              {
 
// #pragma omp for schedule(dynamic, 10)
// #endif
                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++;
                }
 
// #ifdef       USE_OPENMP
//        } // omp parallel
// #endif
 
                // 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:
enum    { NumberOfInputFields = 1 };
typedef LIST< Eagle::tvector3, LIST<Eagle::PhysicalSpace::Metric> > InputTypes;
typedef Fiber::GridInspector< BundleProperty::Anything > GridInspector;
 
        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)
                {}
        };
 
        RefPtr<State> newState() const override
        {
                return new FieldState();
        }
 
 
        IntegralHeart(const string&name, int p, const RefPtr<VCreationPreferences>&VP)
        : VObject(name, p, VP)
        , 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;
 
        bool update(VRequest&R, double precision) override;
 
// experimental AMR code disabled for now, WB 19.JUL.2010
#if     0
        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
                bool apply(const RefPtr<Fiber::Representation>&CartesianLevelRep,
                                    double time, int Level,
                                    const RefPtr<Fiber::Grid>&GridWithCoarsetRefinementLevel,
                                    const RefPtr<ValuePool>&Context) override
                {
                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
                        bool apply (double t, GridID&id, Grid&g) override
                                {
                                        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();
                }
                */
 
        };
#endif
 
 
        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;
 
        VActionNotifier::Progress viP( "Integrating...", max_steps );
                while( hasfinished == false && step_counter < max_steps && viP.proceed( step_counter ) )
                {
                        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 );
#ifdef VERBOSE
//                DebugGridContent(CurrentGrid, "<<<..Initial Iteration...>>>");
#endif
                                GridOperator.storeGrid( grid_name, CurrentGrid, time );
                        }
#ifdef VERBOSE
                        puts("maybe create a copy for the new grid");
#endif
                        NewGrid = GridOperator.advanceGrid( CurrentGrid );
 
#ifdef VERBOSE
            printf("<<<... After NewGrid=GridOperator.advanceGrid(CurrentGrid)...>>>\n");fflush(stdout);
//            DebugGridContent(CurrentGrid, " --- CurrentGrid");
//          DebugGridContent(NewGrid, " --- NewGrid");
            puts("do one integration step of the whole front");
#endif
                        integration_width = GridOperator.getIntegrationSize(CurrentGrid);
                        hasfinished = hasfinished ||
                                      !IntT.advance( CurrentGrid, NewGrid, time, start_count, integration_width );
#ifdef VERBOSE
//           DebugGridContent(NewGrid, "<<<...IntT.advance(CurrentGrid, NewGrid, time, start_count, integration_width)...>>>");
                        puts("refine the grid");
#endif
#ifdef PRT_TIME_ONLY
            printf("<<<<<<<<<<<< timestep = %.8f >>>>>>>>>>>>\n", time); fflush(stdout);
#endif
                        NewGrid = GridOperator.refineGrid( NewGrid );
 
#ifdef VERBOSE
                        puts("fill local data in front");
#endif
                        hasfinished = hasfinished ||
                                      !IntT.extractLocalGridData( NewGrid, time );
#ifdef VERBOSE
//            DebugGridContent(NewGrid, "<<<...After new positions/velocities interpolations...>>>");
                        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();
        }
 
        RefPtr< GridOperatorData<FieldType, LineType> > getGridOperatorData( VRequest&Context )
        {
                return NullPtr();
        }
 
};
 
template<typename FieldType, typename LineType>
bool IntegralHeart<FieldType,LineType>::update(VRequest&Context, double precision)
{
double sec = 0.0;
Timer 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;
 
// extract the emitter grid from field or gridselector, getEmitterGrid is a function defined in IntegralFace.cpp
EmitterFieldData EFD( EmitterGS, EmitterFS, time );
RefPtr<Grid> EmitterGrid = getEmitterGrid( 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.
//
// create a gridname dependent on the grid of the data field and the own module name
string grid_name = FieldSelection.Gridname() + "_" + Name();
 
GridSelector GS;
        myIntegralLines << Context >> GS;
 
// check if a recomputation is necessary, needUpdate is a function defined in IntegralFace.cpp
UpdateCheckData UCD( EmitterGrid, 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);
 
// in case no emitter field or grid is connected, create initial seed points along the diagonal
// of the bounding box of the data fiel
RefPtr<BoundingBox> ToIntegratefieldBBox = getBoundingBox( *FieldSelection.getGrid() );
        if (!ToIntegratefieldBBox) return setStatusError(Context, "No bounding box available.\n");
 
AutoEmitterData AED( EmitterGrid, ToIntegratefieldBBox, 23, FieldSelection.theSourceBundle);
        setDefaultEmissionPoints(Context, AED); // defined in IntegralFace.cpp
 
//
// 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;
 
Eagle::point3 tree_query_scale;
        m_query_scale << Context >> tree_query_scale;
 
RefPtr<LocalFromWorldPoint>
        LocalPointFinder = new LocalFromWorldPoint(*currentSlice.getSlice(), FieldSelection.getGrid(),
                                                   FieldSelection.getGridname(),
                                                   //FieldSelection.getGrid()->getCartesianPositions(),
                                                   UnimapperScale, 0.1, tvector( tree_query_scale ) );
 
double  step_size;
double  min_line_length = 0.0;
double  max_steps = 0.0;
 
        MinLineLength  << Context >> min_line_length;
        MaxSteps       << Context >> max_steps;
        StepSize       << Context >> step_size;
 
        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 = static_cast<unsigned long>(std::floor(max_steps));
 
//Enum AddData, AddMagnitude;
//      IncludeData << Context >> AddData;
//      IncludeMagnitude << Context >> AddMagnitude;
 
        cout << "size: " << EFD.Size << ", seeds: " << AED.nSeeds << endl;
 
unsigned nEmPoints = (EFD.Size==0) ? AED.nSeeds:EFD.Size; // if efd.size is 0 used autocreated seeds instead
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;
 
RefPtr< GridOperatorData<FieldType, LineType> > myGridSetupData = getGridOperatorData( Context );
 
        if(!myGridSetupData)
        {
                puts("IntegralHeart::update(): no grid operator data created!");
//              return setStatusError(Context, "no grid operator data created");
                // in case of Frontstreamlines this is yet null and should be ok
        }
 
GridOperator<FieldType, LineType> integrationGridOperator( FieldSelection.getBundle(),
                                                           myGridSetupData );
 
// TODO make sure emitter grid is there. overrides current version of emitter-grid-point-creation!
RefPtr<Grid> IntegrationGeometry = integrationGridOperator.getPreparedGrid( AED.EmitterGrid,
                                                                            static_cast<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;
 
Options options;
        SomeOptions << Context >> options;
 
        puts("2");
 
        ________________________________________________getTime(sec , Tim_update );
 
        if( Interpol("linear") )
        {
 
        AdvancerLinear  MyAdvancerLinear( FieldSelection, LocalPointFinder, ToIntegratefieldBBox, flowIterator,
                                          step_size, options("dop853"), nPoints, min_line_length, time );
 
                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, options("dop853"), nPoints, min_line_length, time );
 
                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, options("dop853"), nPoints, min_line_length, time );
 
                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::StatusInfo( "Integration computation complete");
 
 
        return true;
}
 
}