IntegralLine.hpp

#include "IntegralFace.hpp"
 
#warning "Including a header from a plugin"
#include <BaseSpace/FieldInterpolator.hpp>
 
#include <eagle/PhysicalSpace.hpp>
//#include <ocean/eagle/STA.hpp>
#include <ocean/shrimp/VEnum.hpp>
#include <ocean/shrimp/VFlagList.hpp>
#include <ocean/shrimp/PhysicalSpace.hpp>
 
#include <memcore/Timer.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 <memcore/Timer.hpp>
 
#include <grid/types/LineSet.hpp>
 
#include <fish/fiber/baseop/RectilinearInterpolation.hpp>
 
#ifdef  USE_OPENMP
#include <omp.h>
#endif
 
/*
  This header file is part of the old streamline/geodesic code:
  and is used by FlowTassel3.cpp and Geodesics.cpp
 
  It should be replaced by IntegralHeart.hpp and FrontStreamlines.cpp,
  FrontPathlines.cpp and FrontGeodesics.cpp.
 
  Since the Geodesics integration was not yet ported to the new structure
  this is still here.
*/
 
 
//#define TIMEIT
//#define VERBOSE
 
namespace FieldLines
{
 
namespace
{
 
// using namespace Wizt;
// using namespace Fiber;
// using namespace Eagle;
// using namespace Eagle::PhysicalSpace;
// using namespace Traum;
 
//typedef       std::string FragmentID_t;
 
typedef RefPtr<FragmentID> 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;
        RefPtr<FragmentID>      fragment;
 
        bool    is_end = false;
 
        tvector      direction;
 
        double  length = 0.0;
 
        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;
        bool inited;
        double scale;
 
        AtomicDataBase( const RefPtr<LocalFromWorldPoint>&LocalPointFinderP,
                        const FieldSelector&FieldSelectionP, const RefPtr<BoundingBox>&BBP,
                        double step_sizeP, double scaleP )//, bool normalizeInversionP, bool doInversionP ) //, const unsigned interpolP )
        : LocalPointFinder(LocalPointFinderP)
        , FieldSelection(FieldSelectionP) // somehow this was not sufficient ?!? see below
        , BB(BBP)
        , step_size(step_sizeP)
        , inited(false)
        , scale( scaleP )
          //, normalizeInversion(normalizeInversionP)
          //, doInversion(doInversionP)
          //, 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!");
        }
 
};
*/
 
 
 
template<typename FieldType, typename LineType, int InterpolType>
class AtomicIntegrator;
 
struct AtomicIntegratorParameterBase
{
        RefPtr<LocalFromWorldPoint> LocalPointFinder;
        FieldSelector FieldSelection;
        RefPtr<BoundingBox> BBP;
        double step_size;
        double scale;
//      bool invNorm, invEnable;
 };
 
template<typename FieldType, typename LineType>
class AtomicIntegratorParameter: public AtomicIntegratorParameterBase
{
};
 
 
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 ]);
#ifdef VERBOSE
                        cout << LineResult[ VertexNr ] << " ";
#endif
                        VertexNr++;
                }
#ifdef VERBOSE
                        cout << endl;
#endif
        }
        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) );
#ifdef VERBOSE
        cout << "CollectLineData(), before setPositions" << endl;
#endif
 
        if( !CartesianVertices.setPositions( LinePositions ) )
          cout << "CollectLineData() ERROR setting Positions" << endl;
 
}
 
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) );
}
 
 
 
// struct       ExtractPointLocationMetric33
// {
//      typedef FieldLines::IntegrationPoint<Eagle::metric33> input_t;
//      typedef Eagle::point3 output_t;
 
//      output_t operator()(const input_t&in) const
//      {
//              return output_t(in.location[0],in.location[1],in.location[2]);
//      }
// };
 
struct  CopyExtractorMetric33
{
        typedef FieldLines::IntegrationPoint<Eagle::metric33> 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::metric33> > > >&Lines, size_t nPoints)
{
RefPtr<Field>   LinePositions = new Field();
        LinePositions->setPersistentData( CollectLines( CopyExtractorMetric33(), 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)
{
#ifdef VERBOSE
  cout << "createLineSet() Lines has lines: " <<  Lines.size() << endl;
#endif
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();
 
#ifdef VERBOSE
        cout << "createLineSet() Line has length: " << LineLength << endl;
#endif
 
                E.resize( LineLength );
                for(index_t i=0; i<LineLength; i++)
                {
                        E[i] = VertexNumber;
#ifdef VERBOSE
                        cout << E[i] << " ";
#endif
                        VertexNumber ++;
                }
#ifdef VERBOSE
                        cout << endl;
#endif
 
        }
 
        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:
//      VActionNotifier::Progress Info;
 
        LineAdvancer(const std::vector<IntegrationResult_t>&LinesP, const bool dopP, const size_t nPointsP,
                     const double min_line_lengthP, const double start_timeP, const AtomicIntegratorParameter<FieldType, LineType> &AIParam)
        : AI(AIParam)
        , Lines(LinesP)
        , dop(dopP)
        , nPoints(nPointsP)
        , min_line_length(min_line_lengthP)
        , start_time(start_timeP)
        , time(start_timeP)
//      , Info("Computing integral line ", Lines.size() )
        {
//                cout << "--------LineAdvancer()" << endl;
//                cout.flush();
 
        }
 
//      LineIntegrator::LineIntegrator(){}
 
 
        size_t getNPoints() { return nPoints;}
 
//      std::vector<IntegrationResult_t>& getLines(){ return Lines; }
 
 
        bool    advanceBreadthFirst( )
        {
#ifdef VERBOSE
                puts("advanceBreadthFirst( ) enter");fflush(stdout);
#endif
 
//              MEMCORE_PROFILE_THIS("LineAdvancer::advanceBreadthFirst( )", 0);
 
                if( dop && !AI.inited )
                        AI.initDop853( Lines.size() );
 
 
        unsigned end_counter = 0;
 
#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(static, 10)
#endif
 
 
                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())) );
//                      VActionNotifier::ProgressInfo("Computing integral line ", p, Lines.size() );
 
//                      if (!Info.proceed(p) )
//                              break;
 
                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;
 
//                      {
//              MEMCORE_PROFILE_THIS("LineAdvancer::advanceBreadthFirst( )", 1);
 
                        if( !dop )
                                test = AI.doEuler( *IntegrationLine, time);
                        else
                                test = AI.doDop853( *IntegrationLine, time, p );
//                      }
 
#ifdef  USE_OPENMP
                        #pragma omp critical
                        {
#endif
                        if( test == false )
                                end_counter++;
                        else
                                nPoints++;
#ifdef  USE_OPENMP
                        } // omp critical
#endif
 
                }
#ifdef  USE_OPENMP
                } // omp parallel
#endif
 
 
#ifdef VERBOSE
                puts("advanceBreadthFirst( ) eexit");fflush(stdout);
#endif
//              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()
        {
//              MEMCORE_PROFILE_THIS("LineAdvancer::finalize()", 0);
 
                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 IntegralLinesParameters
{
  public:
        IntegralLinesParameters(VObject*)
        {}
 
        AtomicIntegratorParameter<FieldType, LineType> getAtomicIntegratorParameter(VRequest&Context)
        { return AtomicIntegratorParameter<FieldType, LineType>(); }
 
};
 
template <typename FieldType, typename LineType>
class   IntegralLines : public IntegralFace, public IntegralLinesParameters<FieldType, LineType>
{
public:
enum    { NumberOfInputFields = 1 };
typedef META::LIST< tvector > 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, tvector >          VecsArray1D_t;
typedef MemArray< 1, RefPtr<FragmentID> >       FragmentIDs_t;
typedef MemArray< 1, BoundingBox >              FragmentBounds_t;
 
 
        TypedSlot<Options>      SpecialOptions;
 
        IntegralLines(const string&name, int p, const RefPtr<VCreationPreferences>&VP)
        : VObject(name, p, VP)
        , Fish<Field>("inputfield")
        , IntegralFace(name, p, VP)
        , IntegralLinesParameters<FieldType, LineType> (this)
#define FORCE_CONSTANT_COORDS_OPTION "ForceConstantCoords"
        , SpecialOptions(this, "options", Options(FORCE_CONSTANT_COORDS_OPTION),2 )
        {}
 
        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;
 
        bool update(VRequest&R, double precision) override;
 
        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;
 
        VActionNotifier::Progress viP( "Computing Integration Line...", max_steps);
                while( hasfinished == false && step_counter < max_steps )
                {
                        if (!viP.proceed(step_counter) )
                        {
                                break;
                        }
                        printProgress(old_nPoints, nPoints, nPointsEst, Context );
                        hasfinished = IntT.advanceBreadthFirst();//,                   //'returns'
                        step_counter++;
                }
                IntT.finalize();//,
//              VActionNotifier::StatusInfo( "Computing Integration Line done.");
                return IntT.getNPoints();
        }
 
};
 
template<typename FieldType, typename LineType>
bool    IntegralLines<FieldType,LineType>::update(VRequest&Context, double precision)
{
//double sec = 0.0;
//VTimer Tim_update;
 
//MEMCORE_PROFILE_THIS("IntegralLines<FieldType,LineType>::update(VRequest&Context, double precision)",0)
 
#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, time );
 
        RefPtr<Field> EmitterCoordinates = getEmitterField( Context, EFD );
 
 
Options options;
        SomeOptions << Context >> options;
 
//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 (options("dop853") ) 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");
 
 
AutoEmitterData AED( EmitterCoordinates, ToIntegratefieldBBox, 23 );
        setDefaultEmissionPoints( Context, AED );
        if (!AED.EmitterCoordinates)
                return setStatusError( Context, "No Emitter Coordinates");
 
 
        //
        // 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 = AED.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;
 
 
Info<Grid> FirstGrid = FieldSelection.BundleSource()->findFirst( FieldSelection.getGridname() );
 
Options mySpecialOptions;
        SpecialOptions << Context >> mySpecialOptions;
 
RefPtr<LocalFromWorldPoint> LocalPointFinder;
 
Eagle::point3 foo;
        foo = {1,1,1};
        m_query_scale << Context >> foo;
 
tvector v_scale = tvector(foo);
 
        cout << v_scale << endl;
//      cin.ignore();
 
        if ( mySpecialOptions(FORCE_CONSTANT_COORDS_OPTION) )
        {
                LocalPointFinder =
                        new LocalFromWorldPoint( *FirstGrid.getSlice(), FirstGrid.getGrid(),
                                                 FieldSelection.getGridname()
                                                 //, FirstGrid.getGrid()->getCartesianPositions(),
                                                 , UnimapperScale, 0.1, v_scale);
        }
        else
        {
                LocalPointFinder =
                        new LocalFromWorldPoint(*currentSlice.getSlice(), FieldSelection.getGrid(),
                                                FieldSelection.getGridname()
                                                //, FieldSelection.getGrid()->getCartesianPositions(),
                                                , UnimapperScale, 0.1, v_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;
 
//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 computation X");fflush(stdout);
//      puts("Marcel was here");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;
}
 
//Enum invNorm;
//Enum invEnable;
 
//      InverseNormalize << Context >> invNorm;
//      InverseEnable << Context >> invEnable;
 
#ifdef VERBOSE
puts("1");
#endif
 
typedef LineAdvancer<FieldType, LineType, FieldInterpolatorBase::LINEAR   > AdvancerLinear;
typedef LineAdvancer<FieldType, LineType, FieldInterpolatorBase::CUBIC    > AdvancerCubic;
typedef LineAdvancer<FieldType, LineType, FieldInterpolatorBase::RECTILINEAR    > AdvancerRectilinear;
typedef LineAdvancer<FieldType, LineType, FieldInterpolatorBase::ANALYTIC > AdvancerAnalytic;
 
int nPointsEst = nEmPoints * (int(max_steps));
AtomicIntegratorParameter<FieldType, LineType> AIP = this->getAtomicIntegratorParameter(Context);
AIP.LocalPointFinder = LocalPointFinder;
AIP.FieldSelection = FieldSelection;
AIP.BBP = ToIntegratefieldBBox;
AIP.step_size = step_size;
AIP.scale = options("invert")?-1.0:1.0;
 
//AIP.invNorm = invNorm("yes");
//AIP.invEnable = invEnable("yes");
 
cout << "integrallines: Fieldselection1: " << AIP.FieldSelection.getFieldName() << endl;
cout << "integrallines: Fieldselection2: " << FieldSelection.getFieldName() << endl;
 
        // main loop
 
#ifdef VERBOSE
        puts("main loop");
//      std::cout << "inversionflags: " << invNorm("yes") << ", " << invEnable("yes") << std::endl;
#endif
//      ________________________________________________getTime(sec , Tim_update );
 
//      {
//      MEMCORE_PROFILE_THIS("IntegralLines<FieldType,LineType>::update(VRequest&Context, double precision)",1)
 
        if( Interpol("linear") )
        {
                if(RefPtr<RectiLinearArray_t>   RectilinearPCoords = FirstGrid.getGrid()->getCartesianPositions()->getFirstData())
                {
        AdvancerRectilinear MyAdvancerRectilinear( Lines, options("dop853"), nPoints, min_line_length, time, AIP );
 
                nPoints = doMainLoop<AdvancerRectilinear>(MyAdvancerRectilinear, max_steps, nPointsEst, nPoints, Context);
 
                }
                else
                {
        AdvancerLinear  MyAdvancerLinear( Lines, options("dop853"), nPoints, min_line_length, time, AIP );
                nPoints = doMainLoop<AdvancerLinear>(MyAdvancerLinear, max_steps, nPointsEst, nPoints, Context);
                }
        }
        else if( Interpol("cubic") )
        {
        AdvancerCubic  MyAdvancerCubic( Lines, options("dop853"), nPoints, min_line_length, time, AIP );
 
                nPoints = doMainLoop<AdvancerCubic>(MyAdvancerCubic, max_steps, nPointsEst, nPoints, Context);
        }
        else
        {
        AdvancerAnalytic  MyAdvancerAnalytic(  Lines, options("dop853"), nPoints, min_line_length, time, AIP );
 
                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
//
//
#ifdef VERBOSE
        puts("IntegralLines::update() 8. Bake");fflush(stdout);
#endif
        {
//      MEMCORE_PROFILE_THIS("IntegralLines<FieldType,LineType>::update(VRequest&Context, double precision) Bake Grid",1)
 
        //
        // 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);
 
        {
//              MEMCORE_PROFILE_THIS("IntegralLines<FieldType,LineType>::update(VRequest&Context, double precision) CollectLineData",1)
                CollectLineData( *(OutputLineSet.CartesianVertices), Lines, nPoints);
                OutputLineSet.Coords = OutputLineSet.CartesianVertices->getPositions();
                assert( OutputLineSet.Coords && "No Coords found");
        }
 
        //
        // 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;
                typedef MemArray<1, FragmentID_t> FragmentIDs_t;
                Ref<FragmentIDs_t> LinesVertFragID( nPoints );
//              MultiArray<1,std::string > LineVFID = *LinesVertFragID;
                MultiArray<1, FragmentID_t>&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( options("addData") )
                {
                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 (options("addMagnitude"))
                {
                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
        //
#ifdef VERBOSE
                puts("before setup standard");fflush(stdout);
#endif
                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"));
 
//      VActionNotifier::StatusInfo( "Integration Line Computation done." );
 
        puts("IntegralLine::update() >>>  EXIT");
 
        return true;
}
 
}
}