|
VISH
0.2
|
#include <ocean/plankton/VModules.hpp> #include <ocean/plankton/VCreator.hpp> #include <bone/FishField.hpp> #include <bone/GridObject.hpp> #include <ocean/GLvish/VRenderObject.hpp> #include <ocean/GLvish/BoundingBox.hpp> #include <ocean/eagle/PhysicalSpace.hpp> #include <ocean/shrimp/VEnum.hpp> #include <field/DirectProductArray.hpp> #include <ocean/plankton/VOperatorCache.hpp> #include <ocean/plankton/VTimer.hpp> #include <fish/fiber/vector/Interpolate.hpp> #include <fish/fiber/vector/LinearIpol.hpp> #include <fish/fiber/baseop/UniGridMapper.hpp> #include <fish/fiber/baseop/ExpandBBox.hpp> #include <ocean/eagle/KDTree.hpp> #include <fish/fiber/field/Fragment.hpp> #include <fish/fiber/grid/Chart.hpp> #include <bone/GridObject.hpp> #include "EmitterPoints.hpp" #include "LinesPoints.hpp" //#define TIMEIT using namespace Wizt; using namespace Fiber; using namespace Eagle; void ________________________________________________getTime(double&sec, VTimer&timer) { #ifdef TIMEIT sec = timer.secs(); #endif } void ________________________________________________printTime(const char*text, const double time) { #ifdef TIMEIT printf("----------------------------------------------------------------------------------- %s:\t%.4f\n", text, time); #endif } class ComputeGridStreamLines : public virtual VRenderObject, public virtual Fish<Grid>, public virtual Fish<Skeleton>, public virtual Fish<Field> { // This class allows to connect an additional grid in a different bundle to the ComputeGridLines node // Used to provide the starting points of the streamlines // Note: Bundle also has its own time parameter struct EmitterGrid : virtual Fish<VObject>, virtual Fish<Bundle>, virtual Fish<Grid> { EmitterGrid( VObject*VO, const Fiber::BundlePtr&BP = Fiber::BundlePtr(), const string&s = "startgrid", int prior=1 ) : Fish<VObject>(VO) , Fish<Bundle>(BP, "startbundle") , Fish<Slice>("starttime", "starttimescale") , Fish<Grid>(s, prior) {} }; public: typedef MemArray<3, Eagle::tvector3> VectorArray_t; typedef MemArray<3, Eagle::point3> CoordsArray_t; typedef MemArray<1, Eagle::point3> CoordsArray1D_t; typedef std::vector<std::vector<point> > Lines_t; typedef DirectProductMemArray<Eagle::point3> ProcArray_t; typedef CreativeIterator<Eagle::point3> CoordsIterator_t; typedef CreativeIterator<double > ScalarIterator_t; typedef CreativeIterator<Eagle::tvector3> VectorIterator_t; typedef std::vector<point> Line_t; typedef pair<point, point> line_t; typedef pair<line_t, rgba_float_t> colorLine_t; typedef pair<point, rgba_float_t> colorPoint_t; struct PerBlockInfo : public ReferenceBase<PerBlockInfo> { point MidPoint; double radius; RefPtr<UniGridMapper> UniMap; RefPtr<FragmentID> FragID; MultiIndex<3> Offset; PerBlockInfo() : ReferenceBase<PerBlockInfo>(this) {} }; struct MyState : State { int RefinementLevel; bool check_all; std::vector<int> candidates; //std::vector<MultiIndex<3> > cand_coord; //std::vector<int > cand_coord_cand; // used for debug drawing only std::vector<colorLine_t> lines; std::vector<colorPoint_t> points; RefPtr<Field> Coords; // Iteraror saves, wheter we have curvilinear or procedural coordinates in the fragments here. bool have_curvi_coords; bool have_fragmented_coords; bool have_fragmented_data; // The field that is used to store the GL cache RefPtr<Field> F; // store midpoints. later this should be exchanged by the kdtree typedef std::vector<RefPtr<PerBlockInfo> > Blocks_t; Blocks_t Blocks; RefPtr<KDTree<3, int> > Tree; double max_radius, last_radius; // remember for skipping recompute int LineLength, Debug, UVWType, TreeAcc; double TreeQueryFac, StepSize; MyState():have_curvi_coords(0), have_fragmented_coords(0), have_fragmented_data(0), max_radius(0.0), last_radius(0.0){} bool changed(const int LineLengthP, const int DebugP, const int UVWTypeP, const int TreeAccP, const double TreeQueryFacP, const double StepSizeP) { if( LineLengthP != LineLength ) { puts("LL changed"); return true; } if( DebugP != Debug ) { puts("DB changed"); return true; } if( UVWTypeP != UVWType ) { puts("UT changed"); return true; } if( TreeAccP != TreeAcc ) { puts("TA changed"); return true; } if( TreeQueryFacP != TreeQueryFac ){ puts("TQ changed"); return true; } if( StepSizeP != StepSize ) { puts("SS changed"); return true; } return false; } void setValues(const int LineLengthP, const int DebugP, const int UVWTypeP, const int TreeAccP, const double TreeQueryFacP, const double StepSizeP) { LineLength = LineLengthP; Debug = DebugP; UVWType = UVWTypeP; TreeAcc = TreeAccP; TreeQueryFac = TreeQueryFacP; StepSize = StepSizeP; } }; // iterator only needed to prepare info for a kd tree and insert info into // the kdtree struct It : FieldFragmentIterator { int c; RefPtr<MyState> S; It(const RefPtr<MyState>_S):c(0), S(_S){} bool apply(const RefPtr<FragmentID>&f, const MemCore::StrongPtr<Fiber::CreativeArrayBase>&DC) { printf("entering iterator\n");fflush(stdout); assert(S); if(!S->Tree) S->Tree = new KDTree<3, int>(); assert(S->Tree); RefPtr<PerBlockInfo> VS = new PerBlockInfo(); // store Fragment id in State for random accessability of fragments in the update functino via the state object VS->FragID = f; // get dimensions of fragment MultiIndex<3> fi, li, lix, liy, liz; bool test; test = makeDimensions( li, DC); li = li[0]-1, li[1]-1, li[2]-1; // <---------------------------------------------- ging nit in Optimize if(!test) { printf("It::apply() Error getting fragment dimsensions"); return false; } // get Coordinate offset of Fragment (in case there's the same Positions field used for multiple data Blocks) MultiIndex<3> Offset; if ( RefPtr<FragmentLocation<3> > F3 = interface_cast<FragmentLocation<3> > (f) ) { Offset = F3->Offset; } else { Offset = 0, 0, 0; } VS->Offset = Offset; // prepare Multiindizes for coordinate data access fi = 0,0,0; fi += Offset; lix = li[0], 0, 0; lix += Offset; liy = 0, li[1], 0; liy += Offset; liz = 0, 0, li[2]; liz += Offset; // store MidPoint and radius of bounding sphere of each fragment, this info is going to be used for storing // fragment ids in a kd tree structure. Alternatively, one could use the less tight but existing bounding info of the fragment. // This should definitely be tried out (performance). tvector p000, p010, p110, p100, p001, p011, p111, p101; // if f == NullPtr data is gotten from the unfragmented field // f == NullPtr if there are no fragments RefPtr<ProcArray_t> PCrds; RefPtr<CoordsArray_t> Crds; // try getting fragmented coords (using the frag id of the data field, which in case should always be consistent) PCrds = S->Coords->getData(f); Crds = S->Coords->getData(f); S->have_fragmented_coords = true; if(!f) S->have_fragmented_coords = false; // if no fragmented found try the unfragmented ones if(!PCrds) { PCrds = S->Coords->getData(); if(PCrds) S->have_fragmented_coords = false; } if(!Crds) { Crds = S->Coords->getData(); if(Crds) S->have_fragmented_coords = false; } //cudaws01.cct.lsu.edu if( Crds ) { p000 = tvector((*Crds)[fi]); p010 = tvector((*Crds)[fi + liy]); p110 = tvector((*Crds)[fi + lix + liy]); p100 = tvector((*Crds)[fi + lix]); p001 = tvector((*Crds)[fi + liz]); p011 = tvector((*Crds)[fi + liy + liz]); p111 = tvector((*Crds)[fi + lix + liy + liz]); p101 = tvector((*Crds)[fi + lix + liz]); S->have_curvi_coords = true; } else { if(!PCrds) { puts("no proc and no curvi coords found"); return false; } p000 = tvector((*PCrds)[fi]); p010 = tvector((*PCrds)[fi + liy]); p110 = tvector((*PCrds)[fi + lix + liy]); p100 = tvector((*PCrds)[fi + lix]); p001 = tvector((*PCrds)[fi + liz]); p011 = tvector((*PCrds)[fi + liy + liz]); p111 = tvector((*PCrds)[fi + lix + liy + liz]); p101 = tvector((*PCrds)[fi + lix + liz]); S->have_curvi_coords == false; } //cout << p000 << p010 << p110 << p100 << p001 << p011 << p111 << p101 << endl; VS->MidPoint = point(0.125 * (p000 + p011 + p110 + p100 + p001 + p011 + p111 + p101)); VS->radius = norm( point(p000) - VS->MidPoint ); S->max_radius = max( S->last_radius, VS->radius ); S->last_radius = VS->radius; S->Blocks.push_back( VS ); cout << "ITERATE inserting " << VS->MidPoint << endl; S->Tree->insert(VS->MidPoint, S->Blocks.size()-1 ); //S->Tree->speak(); c++; return true; } } ; static inline double max(const double a, const double b) { return (a>b)?a:b; } static inline bool pointInBB(point& P, RefPtr<BoundingBox> BBox ) { if( BBox->maxcoord()[0] > P[0] && BBox->mincoord()[0] <= P[0] && BBox->maxcoord()[1] > P[1] && BBox->mincoord()[1] <= P[1] && BBox->maxcoord()[2] > P[2] && BBox->mincoord()[2] <= P[2] ) return true; return false; } //norm class BBSelectorBase { protected: RefPtr<MyState>S; public: BBSelectorBase(const RefPtr<MyState>&_S):S(_S) {assert(S);}; virtual ~BBSelectorBase(){} virtual unsigned int getNext() = 0; virtual bool hasNext() = 0; virtual void reset(const int _last, const point& pos) = 0; }; class BBSelectorMemoryAll : public BBSelectorBase { int size, current, last, count_up, count_down, count; public: BBSelectorMemoryAll(const RefPtr<MyState>_S) : BBSelectorBase(_S), current(0), last(0), count_up(1), count_down(0), count(0) { size = S->Blocks.size(); } unsigned int getNext(); bool hasNext(); void reset(const int _last, const point& pos) { current = 0; last = _last; count = 0; count_up = 1; count_down = 0; } }; class BBSelectorMemoryTree : public BBSelectorBase { int last, current, size; point reset_pos; map<double, int>m; map<double, int>::iterator it; public: BBSelectorMemoryTree(const RefPtr<MyState>_S) : BBSelectorBase(_S), reset_pos(point(0,0,0)) { m.clear(); last=0; current=0; size=0; } unsigned int getNext(); bool hasNext(); void reset( const int _last, const point& pos ); }; class CellTracing { protected: RefPtr<CoordsArray_t> CoordsArr; RefPtr<VectorArray_t> VecsArr; MultiIndex<3> iAx, iAy, iAz; public: CellTracing() { iAx = 1, 0, 0; iAy = 0, 1, 0; iAz = 0, 0, 1; } virtual ~CellTracing(){} void setCoords(const RefPtr<CoordsArray_t> CA) { CoordsArr = CA; } void setVecs (const RefPtr<VectorArray_t> VA) { VecsArr = VA; } inline tvector cross(const tvector&t1, const tvector&t2) { tvector res; res = t1[1]*t2[2] - t1[2]*t2[1], -( t1[0]*t2[2] - t1[2]*t2[0] ), t1[0]*t2[1] - t1[1]*t2[0]; return res; } virtual int pointInCell(const point&P, const int debug, std::vector<colorLine_t>& debuglines) = 0; virtual point localUVW (const point&p, const int debug, std::vector<colorLine_t>& debuglines) = 0; }; class HexaHedralCellNewton : public CellTracing { RefPtr<UniGridMapper> UniMap; point uvw; int point_type; public: HexaHedralCellNewton():CellTracing() {} virtual ~HexaHedralCellNewton(){} void setUniMap ( RefPtr<UniGridMapper> unimap ) { UniMap = unimap; } int pointInCell(const point&P, const int debug, std::vector<colorLine_t>& debuglines); point localUVW (const point&p, const int debug, std::vector<colorLine_t>& debuglines); }; class HexaHedralCellXYTracing : public CellTracing { // objects calced in pointInCell() for use in localUVW() std::vector<point> interval_xy; std::vector<MultiIndex<3> > interval_index_xy, interval_index_xy_perm; MultiIndex<3> Z, uZ; public: HexaHedralCellXYTracing() { interval_xy.resize(4); interval_index_xy.resize(4); interval_index_xy_perm.resize(4); } virtual ~HexaHedralCellXYTracing(){} static inline double dist_xy(const point& p1, const point& p2) { return sqrt( (p2[0]-p1[0])*(p2[0]-p1[0]) + (p2[1]-p1[1])*(p2[1]-p1[1]) ); } // intersect line segment and half-line int testIntersect(const point& p1, const point& p2, const point& p, const tvector& d); int intersect(const point& p1, const point& p2, const point& p, const tvector& d, double*lambda, double*mue); //costs // + 7 * 10 int pointInCell(const point&P, const int debug, std::vector<colorLine_t>& debuglines); // method must only be called after successful pointInCell test! point localUVW (const point&CurrentPoint, const int debug, std::vector<colorLine_t>& debuglines); }; override RefPtr<State> newState() const { return new MyState(); } // this equips the Node with a bundle and a grid slot EmitterGrid StartGrid; TypedSlot<int> LineLength; TypedSlot<double> StepSize; TypedSlot<double> SliceLocation; TypedSlot<Enum> ScalingOptions; TypedSlot<Enum> Debug; TypedSlot<Enum> UVWType; TypedSlot<Enum> Shifter; TypedSlot<Enum> TreeAcc; TypedSlot<double> Scale, Thickness, ArrowAngle, RelativeSize, AtanLength, TreeQueryFac; TypedSlot<int> FirstFragment, LastFragment; TypedSlot<int> nX, nY, nZ; VOutput< ::IntegralLines::LinesPoints> LinesOutput; // TypeList_t AcceptedVectorfieldTypes; VOutput<BundlePtr> myBundle; VOutput<Grid> myIntegralLines; ComputeGridStreamLines(const string&name, int p, const RefPtr<VCreationPreferences>&VP) : VRenderObject(name, p, VP) , Fish<VObject>(this) , Fish<Field>("datafield") , StartGrid(this) , LineLength(this, "linelength", 40) , StepSize(this, "stepsize", 0.02) , SliceLocation(this, "slicecoords", 0.0, 3) , ScalingOptions(this,"normalize", Enum("linear", "norm", "atan"), 3) , Debug(this, "debug", Enum("off", "on")) , UVWType(this, "uvw_type", Enum("XY", "Newton")) , Shifter(this, "shifter", Enum("off", "on")) , TreeAcc(this, "kdtree", Enum("off", "on")) , Scale (this, "scale" , 0.0, 3) , Thickness (this, "thickness" , 1.0 ,NullPtr() , 2) , ArrowAngle (this, "arrowangle" , 45.0 ,NullPtr() , 3) , RelativeSize (this, "headsize" , 0.3 ,NullPtr() , 3) , AtanLength (this, "ascale" , 1.0 ,NullPtr(), 3) , TreeQueryFac(this, "treefac", 1.0) , FirstFragment(this, "firstfragment", 0) , LastFragment (this, "lastfragment" , 10, 4) , nX(this, "iX", 0, 2) , nY(this, "iY", 0, 2) , nZ(this, "iZ", 0, 2) , LinesOutput(self(), "line_output", ::IntegralLines::LinesPoints() ) , myBundle(self(), "outspacetime", BundlePtr() ) , myIntegralLines(self(), "integrallines", GridSelector() ) { printf("ComputeGridStreamLines: I'm Here!!\n");fflush(stdout); FirstFragment.setProperty("max", 2087); LastFragment.setProperty("max", 4000); LineLength.setProperty("max", 500); Scale.setProperty( "min",-10.0); Scale.setProperty( "max", 10.0); Thickness .setProperty( "max", 10.0); ArrowAngle.setProperty( "max", 90.0); AtanLength.setProperty( "min",-10.0); AtanLength.setProperty( "max", 10.0); // AcceptedVectorfieldTypes.insert( FiberType<Eagle::tvector3>::getFiberType().self() ); acceptType<Eagle::tvector3> (); } override bool update(VRequest&R, double precision); override void render(VRenderContext&Context) const; }; //#define VERBOSE bool ComputeGridStreamLines::update(VRequest&Context, double precision) { puts("update");fflush(stdout); VTimer Tim_update; double secs; printf("---- ComputeGridStreamLines::update\n");fflush(stdout); if (!updateSelection(Context) ) { puts("ComputeGridStreamLines::update(): No field!"); removeState(Context); return true; } RefPtr<MyState> S = myState(Context); int ML = getMaxLevel( Context ); BBoxExpander BBE( getBoundingBall( Context ) ); for(int rl = 0; rl <= ML; rl++) { RefPtr<Skeleton> Level = getRefinementLevel(SpacetimeSlot(), Context, rl).second; if (!Level) { printf("ComputeGridStreamLines: update(): No level found!\n"); break; } RefPtr<Representation> LevelRep = Level->DefaultRepresentation(); RefPtr<Field> Coords = LevelRep->Positions(); if (!Coords) { puts("No coord field!"); return true; } S->Coords = Coords; BBE.expand( Coords ); RefPtr<Field> F = getField( LevelRep, Context ); if (!F) { printf("No field for [%s]!\n", Fieldname( Context ).c_str() ); return true; } S->F = F; if(S->F->nFragments() > 1) S->have_fragmented_data = true; else S->have_fragmented_data = false; S->RefinementLevel = rl; printf("We have %d Fragments here!\n", S->F->nFragments() ); // skip line if already done before, iterate over all Blocks and build KDTree if(S->Blocks.size() == 0) { bool test; It MyIt(S); puts("before iterate");fflush(stdout); ________________________________________________getTime(secs, Tim_update); test = F->iterate(MyIt); ________________________________________________printTime("Coords iteration", Tim_update.secs() - secs); if(!test) { puts("iteration failure, abording update"); return false; } printf("xxxxxxxxxxxxxxxxxxxxxxxxx heAr the tree speAK xxxxxxxxxxxxxxxxxxxxxxxxxxxxx\n");fflush(stdout); S->Tree->speak(); } } //S->adjustSliceRange(SliceLocation); setBoundingBall( Context, BBE.BBox); // use caching field for caching the lines pair<double, RefPtr<Fiber::Grid> > inGrid = StartGrid.getGrid(Context); if(!(inGrid.second)) { puts("No emitter grid found"); return true; } // get the infield "Positions" in the a MultiArray, if available, else copy pos into as element into the Array RefPtr< Field > inField; RefPtr<CoordsArray1D_t > EmitterVerts; inField = (*inGrid.second)("Positions"); assert(inField); EmitterVerts = inField->getData(); MultiArray<1, point>&EmitterCrds = *EmitterVerts; puts("got inGrid");fflush(stdout); // Now get the outgrid, if it is available. It can be exctracted from the VOutput // it is used to compare it to the age of the inGrid, to decide if a recalculation // should be done. RefPtr<Grid> outGrid; RefPtr<Field> outField; BundlePtr outBPtr; myBundle << Context >> outBPtr; puts("outGrid1");fflush(stdout); // get the current time of the bundle of the vector field pair< double, RefPtr< Slice > > currentSlice = getSlice(SpacetimeSlot(), Context); printf("Time: %f\n", currentSlice.first ); //get the grid from the VOutput. Bundle pointer and Gridselector have to be used for that purpose GridSelector outGridSelector; myIntegralLines << Context >> outGridSelector; if( !outGridSelector ) { printf("no out Grid\n"); } else outGrid = (*outBPtr)( currentSlice.first , string( (*outGridSelector)() ) ); puts("outGrid2");fflush(stdout); if(!!outGrid) { outField = (*outGrid)("Positions"); assert(outField); } bool valid_outField = !(!outField); printf("got outGrid? Field: %d\n", valid_outField); double step_size, tree_fac; Enum Deb, UV, Shift, TreeA; int line_length = 0; LineLength << Context >> line_length; StepSize << Context >> step_size; Debug << Context >> Deb; UVWType << Context >> UV; Shifter << Context >> Shift; TreeAcc << Context >> TreeA; TreeQueryFac<< Context >> tree_fac; printf("valid_field %d\n", valid_outField); if ( valid_outField ) { printf("in the if %d %d\n", !(S->changed( line_length, Deb(), UV(), TreeA(), tree_fac, step_size )), outField->isNewerThan( *inField )); if ( !(S->changed( line_length, Deb(), UV(), TreeA(), tree_fac, step_size )) && outField->isNewerThan( *inField ) ) { printf("ComputeStreamLines::update(): Line Calc Skipped, Grid is uptodate ...\n"); return true; } } S->setValues(line_length, Deb(), UV(), TreeA(), tree_fac, step_size); puts("2");fflush(stdout); // OCLines->resize( StartPoints.Vertices.size() ); // resize std::vector to number of lines // GLCache::clear( *OCLines ); // if no emitter points available -> quit //if( StartPoints.Vertices.size() <= 0 ) return true; //printf("ComputeGridStreamLines:update(): StartPoint[0]: %f %f %f\n", StartPoints.Vertices[0][0], StartPoints.Vertices[0][1], StartPoints.Vertices[0][2]); S->lines.clear(); //these are just lines drawn for debugging puts("3");fflush(stdout); HexaHedralCellXYTracing HexaXY; HexaHedralCellNewton HexaNewton; BBSelectorMemoryAll BBSelectAll( S ); BBSelectorMemoryTree BBTreeSel( S ); BBSelectorBase* BBSelect; if(TreeA() == 0) BBSelect = &BBSelectAll; else BBSelect = &BBTreeSel; puts("4");fflush(stdout); int last_cand = 0; RefPtr<VectorArray_t> VecsArr; double bound_secs = 0.0; double tree_secs = 0.0; double alllines_secs; unsigned int nEmPoints = EmitterCrds.Size()[0]; ________________________________________________getTime(alllines_secs, Tim_update); unsigned int nth = nEmPoints/10; nth = (nth == 0)?5:nth; puts("5");fflush(stdout); Lines_t compLines; compLines.resize(nEmPoints); for( unsigned int vs = 0; vs < nEmPoints; vs++ ) { unsigned percent; char buf[16]; if( vs % nth == 0) { percent = (unsigned)ceil( (double)vs / nEmPoints * 100); sprintf(buf, "%d %% done", percent); VActionNotifier::ProgressInfo (string(buf)); } printf("ComputeGridStreamLines:update(): Entering StartPointvertex: %d ---------------------------------------------- \n", vs); fflush(stdout); S->candidates.clear(); MultiIndex<1> mvs; mvs = vs; point CurrentPoint = EmitterCrds[mvs]; printf("vs %d, mvs %d, CurrentPOint: %f %f %f\n", vs, mvs[0], CurrentPoint[0], CurrentPoint[1], CurrentPoint[2]); //myLine.push_back(CurrentPoint); compLines[vs].push_back(CurrentPoint); bool abreis_shifted = false; for(int len = 0; len < line_length; len++) { ________________________________________________getTime(secs, Tim_update); int real_candidate = 0; int real_cand_index = 0; int cell_found = 0; //reset block selector double tree_tmp_secs, tmp_secs; ________________________________________________getTime(tree_tmp_secs, Tim_update); BBSelect->reset( last_cand, CurrentPoint ); ________________________________________________getTime(tmp_secs, Tim_update); tree_secs += (tmp_secs - tree_tmp_secs); // printf("Has Next %d\n", BBSelect->hasNext() ); while( BBSelect->hasNext() ) { int j = BBSelect->getNext(); real_cand_index = j; //printf("getNext() %d\n", j);fflush(stdout); //assert(S->Blocks[j]->FragID); // get the bounding box of the fragment RefPtr< CreativeArrayBase > MetaData = S->F->getCreator( S->Blocks[j]->FragID ); assert( MetaData ); double bound_tmp_secs; ________________________________________________getTime(bound_tmp_secs, Tim_update); RefPtr<BoundingBox> BBox; BBox= getFragmentBBox(S->Blocks[j]->FragID, MetaData, *(S->Coords) ); assert( BBox ); ________________________________________________getTime(tmp_secs, Tim_update); bound_secs += (tmp_secs - bound_tmp_secs); VecsArr = S->F->getCreator( S->Blocks[j]->FragID )->create(); assert( VecsArr ); //printf("status of state: have_curvi=%d, have_fragmented_coords=%d, have_fragmented_data=%d\n, ", S->have_curvi_coords, S->have_fragmented_coords, S->have_fragmented_data); // check if point is inside the selected BBox if (BBox->inside(CurrentPoint)) { S->candidates.push_back(j); //printf("ComputeGridStreamLines:update(): index: %d\n", j);fflush(stdout); if( S->have_curvi_coords ) { RefPtr<CoordsArray_t> CoordsArr = S->Coords->getData(S->Blocks[j]->FragID); assert( CoordsArr ); // prepare and invoke search for the cell in the selected block if( UV() == 0 ) { HexaXY.setCoords( ( CoordsArr ) ); HexaXY.setVecs ( ( VecsArr) ); cell_found = HexaXY.pointInCell(CurrentPoint, Deb(), S->lines); } else if( UV() ==1 ) { HexaNewton.setCoords( (CoordsArr) ); HexaNewton.setVecs ( (VecsArr) ); double uni_secs; ________________________________________________getTime(uni_secs, Tim_update); if( !(S->Blocks[ real_cand_index ]->UniMap) ) { S->Blocks[ real_cand_index ]->UniMap = new UniGridMapper( (CoordsArr) ); } ________________________________________________printTime("Unigridmapper:", Tim_update.secs() - uni_secs); HexaNewton.setUniMap( (S->Blocks[ real_cand_index ]->UniMap) ); cell_found = HexaNewton.pointInCell(CurrentPoint, Deb(), S->lines); } } else // do the search for a procedural coords { // direct product arrays are alsways aligned to world xyz coordinates, thus // a positive check of the bounding box guaratnies to be in the coordinate range // of the block. cell_found = true; } last_cand = j; if( cell_found ) { break; } } } ________________________________________________printTime("Find real cand", Tim_update.secs() - secs); ________________________________________________getTime(secs, Tim_update); // printf("cell_found: %d in %d\n", cell_found, real_cand_index); // compute UVW coords and new line point, if a cell was found if( cell_found ) { point float_index; abreis_shifted = false; if(S->have_curvi_coords) { if(UV() == 0) float_index = HexaXY.localUVW(CurrentPoint, Deb(), S->lines); else if(UV() == 1) float_index = HexaNewton.localUVW(CurrentPoint, Deb(), S->lines); } else { RefPtr<ProcArray_t> PCrd; if( S->have_fragmented_coords ) PCrd = S->Coords->getData( S->Blocks[last_cand]->FragID ); else PCrd = S->Coords->getData(); assert(PCrd); float_index = CurrentPoint; bool test = PCrd->locate(float_index); assert(test); // if unfragged coords correct float index by Offset of fragment if( !S->have_fragmented_coords && S->have_fragmented_data) for(unsigned i = 0; i < 3; i++) float_index[i] -= S->Blocks[last_cand]->Offset[i]; } Interpolate<3, tvector, LinearIpol <tvector> > VecField( *(VecsArr), float_index); tvector DataDir = VecField.eval(); //std::cout << "current point: " << CurrentPoint << std::endl; //std::cout << "direction:" << DataDir << std::endl; point OldPoint = CurrentPoint; CurrentPoint = CurrentPoint + step_size * DataDir; if(Deb() == 1) { S->lines.push_back( colorLine_t( line_t( OldPoint, CurrentPoint ), RED) ); } real_candidate = 0; //myLine.push_back(CurrentPoint); compLines[vs].push_back(CurrentPoint); printf("ComputeGridStreamLines:update(): Proceeding to next line point %d in line %d\n\n", len, vs); } else { // Step a half step back if no point was found to ensure, stepping over tight block borders. // This does not seem neccessary at all after correction of an index bug in the UniGridMapper... if( UV() == 1 && abreis_shifted == false && Shift() == 1) { // abreis_shifted = true; // printf("Trying little shift forward instead!\n"); // point a, b; // unsigned max = myLine.size(); // b = myLine[max-1]; // a = myLine[max-2]; // tvector dir = b-a; // b = a + dir * 0.5; // std::cout<< a <<" + 1.2 * " << dir << std::endl; // myLine.pop_back(); // myLine.push_back(b); // CurrentPoint = b; // std::cout<< "CerrentPoint " << CurrentPoint << std::endl; real_candidate = 0; } else { last_cand = 0; printf("ComputeGridStreamLines:update(): No real candidate found, Breaking Line\n"); break; } } ________________________________________________printTime("Compute cand", Tim_update.secs() - secs); } //for len } //for v ________________________________________________printTime("BoundingBox access", bound_secs); ________________________________________________printTime("TreeQuery access", tree_secs); ________________________________________________printTime("All lines comp", Tim_update.secs() - alllines_secs); //::IntegralLines::LinesPoints&computed_lines = (*OCLines); Lines_t&computed_lines = compLines; printf("ComputeStreamLines::update(); computed_lines size %d\n", computed_lines.size()); // LinesOutput << Context << computed_lines; //printf("All Lines Done\nDeleting empyt lines(horrible!)\n"); for( Lines_t::iterator it = computed_lines.begin(); it != computed_lines.end(); it++ ) { if( it->size() == 1) { it->push_back((*it)[0]); } if( it->size() == 0) { it->push_back(point(0,0,0)); it->push_back(point(0,0,0)); } } string grid_name = getGridname(Context) + "_ComputeGridStreamLines"; pair< double, RefPtr< Fiber::Slice > > ActualSlice = getSlice (SpacetimeSlot(), Context); // create a new Grid in the actual slice of the bundle //Grid&LineGrid = (*ActualSlice.second)[grid_name]; Grid&LineGrid = ActualSlice.second->newGrid(grid_name); // Create skeleton for the vertices Skeleton&LineVertices = LineGrid.makeVertices(3); // Create a skeleton ID of dimension 1, and index depth 1 Skeleton&LineEdges = LineGrid[ SkeletonID(1,1) ]; RefPtr<Chart> myCartesian = LineGrid.makeChart( typeid(CartesianChart3D) ); Representation&CartesianVertices = LineVertices[ myCartesian ]; Representation&EdgesAsVertices = LineEdges[ LineVertices ]; RefPtr<Field> LineCoords = CartesianVertices[ FIBER_POSITIONS ]; RefPtr<Field> LineSetEdges = EdgesAsVertices [ FIBER_POSITIONS ]; MultiIndex<1> NumberOfLineVertices(0); { BundlePtr BPtr; findSpacetimeSlot() << Context >> BPtr; myBundle << Context << BPtr; GridSelector GS( new StringSelection( grid_name ) ); myIntegralLines << Context << GS; } ________________________________________________getTime(secs, Tim_update); for( index_t ii = 0; ii < computed_lines.size(); ii++ ) { MultiIndex<1> computed_linesInd( computed_lines[ii].size() ); NumberOfLineVertices += computed_linesInd; printf("line %d, has %d verteces\n", ii, computed_lines[ii].size()); } printf( "Number of Vertices: %d\n", NumberOfLineVertices[0] ); fflush(stdout); typedef MemArray<1, point> LinesCoordsArr_t; RefPtr<LinesCoordsArr_t> LinesCoordsData = new LinesCoordsArr_t(NumberOfLineVertices); MultiArray<1,point > LineCrds = *LinesCoordsData; { index_t vert_count = 0; for(index_t lin = 0; lin < computed_lines.size(); lin++ ) { printf("\nEntering line %d of %d\n", lin, computed_lines.size());fflush(stdout); for(index_t ver = 0; ver < computed_lines[lin].size() ; ver++ ) { printf(" V(%d %d %d) ", ver, computed_lines[lin].size(), lin);fflush(stdout); printf(" C(%f %f %f)\n", computed_lines[lin][ver][0], computed_lines[lin][ver][1], computed_lines[lin][ver][2]); MultiIndex<1> n( vert_count ); LineCrds[ n ] = computed_lines[lin][ver]; vert_count++; } } //printf("\n"); } LineCoords->setPersistentData( LinesCoordsData ); // printf("Coords done!\n"); typedef MemArray<1, std::vector<index_t> > EdgesArray_t; RefPtr<EdgesArray_t> EdgesArray = new EdgesArray_t( computed_lines.size() ); MultiArray<1, std::vector<index_t> >&LinesEdges = *EdgesArray; { index_t vert_count = 0; for(index_t lin = 0; lin < computed_lines.size(); lin++ ) { //printf("\nEntering line %d of %d\n", lin, computed_lines.Lines.size());fflush(stdout); std::vector<index_t>&E = LinesEdges[ lin ]; E.resize( computed_lines[lin].size() ); //printf("2\n");fflush(stdout); for(index_t ver = 0; ver < computed_lines[lin].size() ; ver++ ) { //printf("(E(%d %d)\n", ver, vert_count);fflush(stdout); E[ver] = vert_count; vert_count++; } } } LineSetEdges->setPersistentData( EdgesArray ); ________________________________________________printTime("Copy in Bundle", Tim_update.secs() - secs); // the bundle data object has to be touched, to be updated correctly if(RefPtr<VObject> vobj = SpacetimeSlot()->Source() ) vobj->touch(); ________________________________________________printTime("update Time", Tim_update.secs()); printf("ComputMultipleStreamLines::update() Written into Bundle, setPersitentData and touched Spacetimeslot()\n"); VActionNotifier::ProgressInfo (string("Streamline computation complete")); return true; } static inline void glLine(const ComputeGridStreamLines::point&A, const ComputeGridStreamLines::tvector&dir) { // glNormal( *dir ); glVertex( A ); glVertex( A+dir ); } void ComputeGridStreamLines::render(VRenderContext&Context) const { /* Enum Deb; Debug << Context >> Deb; if(Deb() == 1) { RefPtr<FieldState> state = getState(Context); if (!state) return; if (state->Blocks.size() < 1) { puts("override void ComputeGridStreamLines::render(VRenderContext&Context) const: NOTHING to render..."); return; } assert(state); int count = 0; glBegin(GL_POINTS); for( std::vector<colorPoint_t>::iterator it = state->points.begin(); it != state->points.end(); it++) { glColor3f((it->second)[0], (it->second)[1], (it->second)[2]); glVertex(it->first); count++; } glEnd(); //printf("%d points drawed\n", count); glBegin(GL_LINES); for( std::vector<colorLine_t>::iterator it = state->lines.begin(); it != state->lines.end(); it++) { glColor3f((it->second)[0], (it->second)[1], (it->second)[2]); glVertex((it->first).first); glVertex((it->first).second); } glEnd(); } // render candidate boxes, only candidates of one line are rendered, of the last one. // So, if the last line is not in data range (only emit point visible) no candidate boxes // are visible. // for(unsigned int i = 0; i < myState->candidates.size(); i++ ) // { // RefPtr<BoundingBox> BBox = myState->Blocks[myState->candidates[i]]->BBox; // glColor3f(0.4,0.4,0.4); // glDrawBox( point(BBox->maxcoord()), point(BBox->mincoord()) ); // } //printf("Render: draw box ? sizes: %d %d\n", myState->cand_coord_cand.size() , myState->cand_coord.size()); // render cand_intervals //printf("ComputeGridStreamLines:render(): sizes: cand_coord %d, cand_coord_cand %d\n", myState->cand_coord.size(), myState->cand_coord_cand.size()); // for(unsigned int j = 0; j < myState->cand_coord_cand.size(); j = j+4 ) // { // MultiIndex<3> m1, m2, m3, m4; // m1 = myState->cand_coord[j]; // m2 = myState->cand_coord[j+1]; // m3 = myState->cand_coord[j+2]; // m4 = myState->cand_coord[j+3]; // //printf("Render: draw box ix %d iy %d in %d\n", m1[0], m1[1], myState->cand_coord_cand[j]); // glBegin( GL_LINES ); // glVertex( (*(myState->Blocks[ myState->cand_coord_cand[j] ]->CoordsArr))[ m1 ] ); // glVertex( (*(myState->Blocks[ myState->cand_coord_cand[j] ]->CoordsArr))[ m2 ] ); // glVertex( (*(myState->Blocks[ myState->cand_coord_cand[j] ]->CoordsArr))[ m2 ] ); // glVertex( (*(myState->Blocks[ myState->cand_coord_cand[j] ]->CoordsArr))[ m3 ] ); // glVertex( (*(myState->Blocks[ myState->cand_coord_cand[j] ]->CoordsArr))[ m3 ] ); // glVertex( (*(myState->Blocks[ myState->cand_coord_cand[j] ]->CoordsArr))[ m4 ] ); // glVertex( (*(myState->Blocks[ myState->cand_coord_cand[j] ]->CoordsArr))[ m4 ] ); // glVertex( (*(myState->Blocks[ myState->cand_coord_cand[j] ]->CoordsArr))[ m1 ] ); // glEnd(); // } // // glLineWidth(2); // glBegin( GL_LINES ); // for(unsigned int j = 0; j < myState->lines.size(); j+=2 ) // { // glColor3f(0, 1, 1); // glVertex ( myState->lines[j]); // glVertex ( myState->lines[j+1]); // } // glEnd(); // glPointSize(4); // glColor3f(1.0, 1.0, 0.0); // glBegin( GL_POINTS ); // for(unsigned int j = 0; j < myState->lines.size(); j+=2 ) // { // glColor3f(0, 1, 1); // glVertex ( myState->lines[j]); // } // glEnd(); // int nx, ny, nz; // nX << Context >> nx; // nY << Context >> ny; // nZ << Context >> nz; // MultiIndex<3> bla; // bla = nx, ny, nz; // glColor3f(0, 1.0, 1.0); // glBegin(GL_POINTS); // glVertex( (*(myState->Blocks[ First ]->CoordsArr))[ bla ] ); // glEnd(); // // glBegin( GL_POINTS ); // // for(unsigned int j = 0; j < myState->points.size(); j++ ) // // { // // glColor3f(1, 1, 0); // // glVertex(myState->points[j]); // // } // // glEnd(); // // puts(" render ok"); // } // glDisable(GL_LINE_SMOOTH); // } */ } int ComputeGridStreamLines::HexaHedralCellXYTracing::testIntersect(const point& p1, const point& p2, const point& p, const tvector& d) { tvector d_seg = p2 - p1; double lambda = 0; double mue = 0; lambda = ( p[1]*d[0] + p1[0] * d[1] - p[0] * d[1] - p1[1]*d[0]) / ( d_seg[1]*d[0] - d_seg[0]*d[1] ); if( lambda < 0.0 || lambda > 1.0 ) return 0; mue = (p1[0] + lambda * d_seg[0] - p[0]) / d[0]; if( mue < 0.0 ) return 0; return 1; } int ComputeGridStreamLines::HexaHedralCellXYTracing::intersect(const point& p1, const point& p2, const point& p, const tvector& d, double*lambda, double*mue) { tvector d_seg = p2 - p1; *lambda = 0; *mue = 0; *lambda = ( p[1]*d[0] + p1[0] * d[1] - p[0] * d[1] - p1[1]*d[0]) / ( d_seg[1]*d[0] - d_seg[0]*d[1] ); if( *lambda < 0.0 || *lambda > 1.0 ) return 0; *mue = (p1[0] + *lambda * d_seg[0] - p[0]) / d[0]; if( *mue < 0.0 ) return 0; return 1; } int ComputeGridStreamLines::HexaHedralCellXYTracing::pointInCell(const point&P, const int debug, std::vector<colorLine_t>& debuglines) { MultiIndex<3> FieldSize = CoordsArr->Size(); //printf("ComputeGridStreamLines:update(): Candidate size: %d %d %d\n", FieldSize[0], FieldSize[1], FieldSize[2] ); MultiIndex<3> ifirst, ilast, ifx, ify, ifz; point float_index, max_float_index, p1, p2; unsigned int ix, iy; max_float_index = FieldSize[0] - 1, FieldSize[1] - 1, FieldSize[2] - 1; ilast = FieldSize[0] - 1, FieldSize[1] - 1, FieldSize[2] - 1; ifirst = 0,0,0; ifx = 1, 0, 0; ify = 0, 1, 0; ifz = 0, 0, 1; p1 = (*CoordsArr)[ifirst]; p2 = (*CoordsArr)[ilast]; tvector vx, vy, vz; vx = (*CoordsArr)[ifx] - p1; vy = (*CoordsArr)[ify] - p1; vz = (*CoordsArr)[ifz] - p1; // if ilast... < 1 --> break // printf("vx: %f %f %f \n", vx[0], vx[1], vx[2]); // printf("vy: %f %f %f \n", vy[0], vy[1], vy[2]); // printf("vz: %f %f %f \n", vz[0], vz[1], vz[2]); // here we are looking for the index that represents the z-Axis, since we are only looking at // the xy coords for the cell finding later. // re-index-index Z and un-re-index-index uZ for reindexing und unreindexing z-index to the 3rd index,of MultiIndizes, // which happens not to be the index for the z-direction generally. if( vx[2] != 0.0) { Z = 2, 1, 0; uZ = 2, 1, 0; //printf("x is z\n"); } else if( vy[2] != 0.0) { Z = 2, 0, 1; uZ = 1, 2, 0; //printf("y is z\n"); } else if( vz[2] != 0.0) { Z = 0, 1, 2; uZ = 0, 1, 2; //printf("z is z\n"); } else { Z = 0, 1, 2; uZ = 0, 1, 2; //printf("nothing is z\n"); } //printf("p1: %f %f %f, p2: %f %f %f\n", p1[0], p1[1], p1[2], p2[0], p2[1], p2[2]); tvector testdir(0.77753, 0.2274, 0.111); //tvector testdir(0.111, 0.45, 0); // only x and y coords are crucial to decide whether the point is really contained (cylinrical coords) // find by intersection of an line in z-plane intersection the bouning rectangle of a coordinate interval in the z-plane //printf("Testing Intervals-------: %d %d %d\n", ilast[0], ilast[1], ilast[2]); //printf("Flipped Intervals-------: %d %d %d\n", ilast[Z[0]], ilast[Z[1]], ilast[Z[2]]); int real_candidate; for(ix = 0; ix < ilast[ Z[0] ]; ix++) { for( iy = 0; iy < ilast[ Z[1] ]; iy++) { //printf("Testing Interval: %d of %d, %d of %d\n", ix, ilast[ Z[0]], iy, ilast[ Z[1] ] ); //fflush(stdout); real_candidate = 0; interval_index_xy[0] = ix, iy, 0; interval_index_xy[1] = ix+1, iy, 0; interval_index_xy[2] = ix+1, iy+1, 0; interval_index_xy[3] = ix, iy+1, 0; //printf("1");fflush(stdout); //for(int i = 0; i < 4; i++) interval_index_xy_perm[0] = interval_index_xy[0][ uZ[0] ], interval_index_xy[0][ uZ[1] ], interval_index_xy[0][ uZ[2] ]; interval_index_xy_perm[1] = interval_index_xy[1][ uZ[0] ], interval_index_xy[1][ uZ[1] ], interval_index_xy[1][ uZ[2] ]; interval_index_xy_perm[2] = interval_index_xy[2][ uZ[0] ], interval_index_xy[2][ uZ[1] ], interval_index_xy[2][ uZ[2] ]; interval_index_xy_perm[3] = interval_index_xy[3][ uZ[0] ], interval_index_xy[3][ uZ[1] ], interval_index_xy[3][ uZ[2] ]; //printf("2");fflush(stdout); //for(int i = 0 ; i < 4 ; i ++) // printf("permindex[%d]: %d %d %d\n", i, interval_index_xy_perm[i][0], interval_index_xy_perm[i][1], interval_index_xy_perm[i][2]); interval_xy[0] = (*CoordsArr)[interval_index_xy_perm[0]]; //t1 interval_xy[1] = (*CoordsArr)[interval_index_xy_perm[1]]; //t2 interval_xy[2] = (*CoordsArr)[interval_index_xy_perm[2]]; //t3 interval_xy[3] = (*CoordsArr)[interval_index_xy_perm[3]]; //t4 //printf("3");fflush(stdout); real_candidate += testIntersect(interval_xy[0], interval_xy[1], P, testdir); real_candidate += testIntersect(interval_xy[1], interval_xy[2], P, testdir); real_candidate += testIntersect(interval_xy[2], interval_xy[3], P, testdir); real_candidate += testIntersect(interval_xy[3], interval_xy[0], P, testdir); //printf("4");fflush(stdout); // S->lines.push_back( interval_xy[0] ); S->lines.push_back( interval_xy[1] ); // S->lines.push_back( interval_xy[1] ); S->lines.push_back( interval_xy[2] ); // S->lines.push_back( interval_xy[2] ); S->lines.push_back( interval_xy[3] ); // S->lines.push_back( interval_xy[3] ); S->lines.push_back( interval_xy[0] ); if(real_candidate == 1) { //printf("REAk CAND found at index %d %d %d\n", interval_index_xy_perm[0][0], interval_index_xy_perm[0][1], interval_index_xy_perm[0][2]); if(debug == 1) { debuglines.push_back( colorLine_t( line_t( interval_xy[0], interval_xy[1]), GREEN) ); debuglines.push_back( colorLine_t( line_t( interval_xy[1], interval_xy[2]), GREEN) ); debuglines.push_back( colorLine_t( line_t( interval_xy[2], interval_xy[3]), GREEN) ); debuglines.push_back( colorLine_t( line_t( interval_xy[3], interval_xy[0]), GREEN) ); } return 1; } } } return 0; } point ComputeGridStreamLines::HexaHedralCellXYTracing::localUVW (const point&CurrentPoint, const int debug, std::vector<colorLine_t>& debuglines) { if( !VecsArr) return point(0,0,0); if( !CoordsArr) return point(0,0,0); if( interval_xy.size() != 4) return point(0,0,0); if( interval_index_xy.size() != 4) return point(0,0,0); if( interval_index_xy_perm.size() != 4) return point(0,0,0); //printf("ComputeGridStreamLines:update(): YEES real cand found at %d ,%d in %d\n", ix, iy, real_cand_index); //printf("ComputeGridStreamLines:update(): sizes: %d %d \n", S->cand_coord.size(), S->cand_coord_cand.size() ); // decide in wich triangle the point is: tvector r = interval_xy[2] - interval_xy[0]; double k = r[1]/r[0]; double d = -k * interval_xy[0][0] + interval_xy[0][1]; std::vector<int> in_eras; // delete via reindexing, to make it a bit confuseing in_eras.resize(4); in_eras[0] = 0; in_eras[1] = 1; in_eras[2] = 2; in_eras[3] = 3; //S->points.push_back(point(CurrentPoint[0], CurrentPoint[1], interval_xy[0][2] )); // is point over line? if( CurrentPoint[1] > k * CurrentPoint[0] + d ) { if( interval_xy[1][1] > k * interval_xy[1][0] + d ) { // interval point 1 is also over -> erase point 3 in_eras.erase(in_eras.begin() + 3); // S->points.push_back( interval_xy[3] ); } else { in_eras.erase(in_eras.begin() + 1); // S->points.push_back( interval_xy[1] ); } } else // is under line { if( interval_xy[1][1] <= k * interval_xy[1][0] + d ) { // interval point 1 is also under -> erase point 3 in_eras.erase(in_eras.begin() + 3); // S->points.push_back( interval_xy[3] ); } else { in_eras.erase(in_eras.begin() + 1); // S->points.push_back( interval_xy[1] ); } } //S->lines.push_back( interval_xy[ in_eras[0] ] ); //S->lines.push_back( interval_xy[ in_eras[1] ] ); //S->lines.push_back( interval_xy[ in_eras[1] ] ); //S->lines.push_back( interval_xy[ in_eras[2] ] ); //S->lines.push_back( interval_xy[ in_eras[2] ] ); //S->lines.push_back( interval_xy[ in_eras[0] ] ); //in_eras.erase(in_eras.begin() + dist_max_index); //printf("ineras: %d %d %d, size %d\n", in_eras[0], in_eras[1], in_eras[2], in_eras.size()); // interpol via cutting lines in z plane point pv; pv = CurrentPoint[0], CurrentPoint[1], interval_xy[ in_eras[0] ][2]; tvector dir = pv - interval_xy[ in_eras[0] ]; double lambda, mue; intersect( interval_xy[ in_eras[1] ], interval_xy[ in_eras[2] ], CurrentPoint, dir, &lambda, &mue ); //printf("lambda %f, mue %f\n", lambda, mue); if(debug == 1) { //S->points.push_back( interval_xy[ in_eras[0] ] + mue * dir ); //S->points.push_back( interval_xy[ in_eras[0] ]); //S->points.push_back( interval_xy[ in_eras[1] ] + lambda * (interval_xy[ in_eras[2] ] - interval_xy[ in_eras[1] ]) ); } tvector tmp_pos= lambda* tvector( interval_index_xy_perm[ in_eras[1] ][0], interval_index_xy_perm[ in_eras[1] ][1], interval_index_xy_perm[ in_eras[1] ][2] ) + (1 - lambda)* tvector( interval_index_xy_perm[ in_eras[2] ][0], interval_index_xy_perm[ in_eras[2] ][1], interval_index_xy_perm[ in_eras[2] ][2] ); point s = interval_xy[ in_eras[1] ] + lambda * (interval_xy[ in_eras[2] ] - interval_xy[ in_eras[1] ]); double n1 = dist_xy( s, pv ); double n2 = dist_xy( pv, interval_xy[ in_eras[0] ] ); double ng = n1 + n2; n1 = n1 / ng; point float_index; float_index = n1 * tmp_pos + (1 - n1)* tvector( interval_index_xy_perm[ in_eras[0] ][0], interval_index_xy_perm[ in_eras[0] ][1], interval_index_xy_perm[ in_eras[0] ][2] ); // printf("lambda %f, mue %f, tmp_pos %f %f\n", lambda, mue, float_index[0], float_index[1]); MultiIndex<3> FieldSize = CoordsArr->Size(); MultiIndex<3> ifirst, ilast; ilast = FieldSize[0] - 1, FieldSize[1] - 1, FieldSize[2] - 1; ifirst = 0,0,0; point p1, p2, max_float_index; p1 = (*CoordsArr)[ifirst]; p2 = (*CoordsArr)[ilast]; max_float_index = ilast[0], ilast[1], ilast[2]; // calc accurate z-float index (cheating by assuming kind of cylindrical coordinates with z as axis) double z_height = p2[2] - p1[2]; // interpolated float index float_index[ Z[2] ] = (CurrentPoint[2] - p1[2]) / z_height * max_float_index[Z[2]] ; //printf("z_height: %f, Current_z - p1z: %f, max_ind: %f\n", z_height, (CurrentPoint[2] - p1[2]), max_float_index[Z[2]]); //printf("ComputeGridStreamLines:update(): index: %f, %f, %f\n", float_index[0], float_index[1], float_index[2] ); MultiIndex<3> c0, c1, c2, c3, c4, c5, c6, c7; c0 = (int)round(floor(float_index[0])), (int)round(floor(float_index[1])), (int)round(floor(float_index[2])); for(int i = 0; i < (int)round(max_float_index[Z[2]]) - 1; i++ ) { MultiIndex<3> da, db; da = c0[0], c0[1], c0[2]; db = c0[0], c0[1], c0[2]; da[Z[2]] = i; db[Z[2]] = i+1; point pa, pb; pa = (*CoordsArr)[ da ]; pb = (*CoordsArr)[ db ]; //printf("da %d %d %d, db %d %d %d\n", da[0], da[1], da[2], db[0], db[1], db[2]); //printf("dists [%d]: norm %f, z-dist: %f\n", i, norm((pb-pa)), pb[2]-pa[2]); //printf("Its %f not here! betweeb %f and %f\n", CurrentPoint[2], pa[2], pb[2]); if( pa[2] <= CurrentPoint[2] && CurrentPoint[2] <= pb[2] ) { //printf("Its heeeeeere! betweeb %f and %f\n", pa[2], pb[2]); break; } if (pa[2] >= CurrentPoint[2] && CurrentPoint[2] >= pb[2] ) { //printf("Its heeeeeere! betweeb %f and %f\n", pa[2], pb[2]); float_index[Z[2]] = 1 - (CurrentPoint[2] - pb[2]) / (pa[2] - pb[2]) + i; //printf("zind: %f\n", float_index[Z[2]]); break; } } //printf("ComputeGridStreamLines:update(): index: %f, %f, %f\n", float_index[0], float_index[1], float_index[2] ); if(debug == 1) { c0[Z[2]] = (int)round(floor(float_index[Z[2]])); //printf("c0: %d %d %d\n", c0[0], c0[1], c0[2]); c1 = c0[0], c0[1]+1, c0[2]; c2 = c0[0]+1, c0[1]+1, c0[2]; c3 = c0[0]+1, c0[1], c0[2]; c4 = c0[0], c0[1], c0[2]+1; c5 = c0[0], c0[1]+1, c0[2]+1; c6 = c0[0]+1, c0[1]+1, c0[2]+1; c7 = c0[0]+1, c0[1], c0[2]+1; debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c0], (*CoordsArr)[c1]), YELLOW) ); debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c1], (*CoordsArr)[c2]), YELLOW) ); debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c2], (*CoordsArr)[c3]), YELLOW) ); debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c3], (*CoordsArr)[c0]), YELLOW) ); debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c4], (*CoordsArr)[c5]), YELLOW) ); debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c5], (*CoordsArr)[c6]), YELLOW) ); debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c6], (*CoordsArr)[c7]), YELLOW) ); debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c7], (*CoordsArr)[c4]), YELLOW) ); debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c0], (*CoordsArr)[c4]), YELLOW) ); debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c1], (*CoordsArr)[c5]), YELLOW) ); debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c2], (*CoordsArr)[c6]), YELLOW) ); debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c3], (*CoordsArr)[c7]), YELLOW) ); } return float_index; } int ComputeGridStreamLines::HexaHedralCellNewton::pointInCell(const point&P, const int debug, std::vector<colorLine_t>& debuglines) { if(!UniMap) { printf("No unimap available");fflush(stdout); return -1; } point_type = UniMap->localCellCoordinatesFromCurviGrid( P, uvw ); return point_type; } point ComputeGridStreamLines::HexaHedralCellNewton::localUVW (const point&p, const int debug, std::vector<colorLine_t>& debuglines) { return uvw; } unsigned int ComputeGridStreamLines::BBSelectorMemoryTree::getNext() { if(current == 0) { // printf("returned last, without tree\n"); current++; return last; } if(current == 1) // query the tree here { //printf("querying creating tree\n"); m.clear(); S->Tree->nearest_range(reset_pos, S->max_radius * S->TreeQueryFac, m); //printf("Map Size %d\n", m.size()); // for(map<double,int>::iterator i = m.begin(); i != m.end(); i++) // printf("map: %f %d\n", i->first ,i->second); size = m.size()+1; current++; if(m.size() == 0) { return 0; } it = m.begin(); int tmp = it->second; it++; return tmp; } else { int tmp = it->second; it++; current++; return tmp; } } void ComputeGridStreamLines::BBSelectorMemoryTree::reset(const int _last, const point& pos) { last = _last; reset_pos = pos; current = 0; it = m.begin(); } bool ComputeGridStreamLines::BBSelectorMemoryTree::hasNext() { if(current == 0) // theres always a last to test return true; if(current == 1) // tree query has to be done return true; if( current < size) return true; else return false; } unsigned int ComputeGridStreamLines::BBSelectorMemoryAll::getNext() { if( count % 2 == 0 ) { if( (last + count_up) < size ) { current = last + count_up; count_up++; } else { current = last - count_down; count_down++; } } else { if( (last - count_down) >= 0 ) { current = last - count_down; count_down++; } else { current = last + count_up; count_up++; } } count++; return current; } bool ComputeGridStreamLines::BBSelectorMemoryAll::hasNext() { if( count < size ) return true; else return false; } static Ref<VCreator<ComputeGridStreamLines, AcceptList<Fiber::Field> > > MyComputeGridStreamLinesCreator("IntegralLines/ComputeGridStreamLines");
1.7.5