LocalFromWorldPoint.hpp

#ifndef __LOCALFROMWORLDPOINT_HPP
#define __LOCALFROMWORLDPOINT_HPP

#include <fish/fiber/field/DirectProductArray.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 <fish/fiber/field/Fragment.hpp>
#include <fish/fiber/grid/Chart.hpp>

#include <fish/fiber/baseop/BundleBaseSpaceFunctions.hpp>

#include <eagle/PhysicalSpace.hpp>
#include <aerie/KDTree.hpp>
#include <eagle/FixedArray.hpp>

#include <memcore/RefPtr.hpp>
#include "FragmentSkeleton.hpp"
#include "PointSearch.hpp"


#include "gridopDllApi.h"

namespace Fiber
{

using namespace Eagle;
using namespace Eagle::PhysicalSpace;


class   gridop_API LocalFromWorldPoint : public MemCore::ReferenceBase<LocalFromWorldPoint >
{
        // info retrieved from iterator, which is called once
//      bool curvi_coords; 
//      bool fragmented_coords; 

//      bool inited; 

//      RefPtr<Slice> slice;
        RefPtr<Field> coords; 

        CartesianFragments      FragmentProperties; 
        PointSearch             UniformPointSearch;

        RefPtr<Grid>&grid() { return FragmentProperties.TheGrid; }


        double TreeQueryFac;
        RefPtr<KDTree<3, int> > MyTree;

        static  inline MultiIndex<3> getFragmentOffset(const RefPtr<FragmentID>&f)
        {
         MultiIndex<3> Offset; 
                 if ( RefPtr<FragmentLocation<3> > F3 = interface_cast<FragmentLocation<3> > (f) )
                 {
                         Offset = F3->Offset; 
                 } 
                 else
                 {
                         Offset = 0, 0, 0;
                 } 

                 return Offset;
        }

        typedef MemArray<3, Eagle::tvector3>           VectorArray_t; 
        typedef MemArray<3, double>                    ScalarArray_t; 
        typedef MemArray<3, Eagle::point3>             CoordsArray_t; 
        typedef MemArray<1, Eagle::point3>             CoordsArray1D_t;
        typedef MemArray<1, std::string>               FragmentIDs_t; 

        typedef FixedArray<point,8>                    Vertexlist_t;

        typedef MemArray<1, RefPtr<BoundingBox> >      FragmentBounds_t;
        typedef DirectProductMemArray<Eagle::point3>   ProcArray_t; 

        // Remember Fragment IDs and Bounding Boxes, and UniGRid object in class after initialization 
//      RefPtr<FragmentIDs_t>    MyFragIds; 
//      RefPtr<FragmentBounds_t> MyFragBBs; 
        RefPtr<Field>&UniMapperField()
        {
                return UniformPointSearch.UniMapperField;
        }

        struct  MyIterator : Field::Iterator
        {
                int             c;

                RefPtr<Field> FragmentBoxes;

                RefPtr<Field>           Coords;
                RefPtr<KDTree<3, int> > Tree;
                double last_radius;

                bool build_tree;

//              RefPtr<FragmentIDs_t>    MyFragIds; 
//              RefPtr<FragmentBounds_t> MyFragBBs;

                MyIterator(const RefPtr<Field>&FragmentBoxesPositions,
                           const RefPtr<Field>&CoordsP, 
                           const RefPtr<KDTree<3,int> >&TreeP,
                           const bool build_treeP)
//, const RefPtr<FragmentIDs_t>&FragIdsP, const RefPtr<FragmentBounds_t>&FragBBsP )
                :c(0)
                , FragmentBoxes( FragmentBoxesPositions )
                , Coords(CoordsP)
                , Tree(TreeP)
                , last_radius(0)
                , build_tree(build_treeP)
//              , MyFragIds(FragIdsP)
//              , MyFragBBs(FragBBsP)
                {}
                
                bool apply(const RefPtr<FragmentID>&f, const MemCore::StrongPtr<Fiber::CreativeArrayBase>&DC);

        } ; 



        class BBSelectorBase
        {
        protected:

                
        public:
                BBSelectorBase() {}
                virtual ~BBSelectorBase(){}
        
                virtual unsigned int getNext() = 0;
                virtual bool hasNext() = 0;
                virtual void reset(const int _last, const point& pos) = 0;
        };      

 public:
        class BBSelectorMemoryTree : public BBSelectorBase
        {
                int last, current, size;
                point reset_pos;
                map<double, int>m;
                map<double, int>::iterator it;
                double TreeQueryFac;
                RefPtr<KDTree<3,int> > Tree;
        public:
                BBSelectorMemoryTree(const RefPtr<KDTree<3, int> >TreeP, const double TreeQueryFacP) 
                : reset_pos(point(0,0,0))
                , TreeQueryFac(TreeQueryFacP)
                , Tree(TreeP)
                {
                        m.clear(); 
                        last=0; 
                        current=0; 
                        size=0;
                }
                
                unsigned int getNext();
                bool hasNext();
                void reset( const int _last, const point& pos );
        };
        
        
        class FindCell
        {
        protected:
                RefPtr<CoordsArray_t>          CoordsArr;
                RefPtr<VectorArray_t>          VecsArr;
                MultiIndex<3> iAx, iAy, iAz; 
        public:
                FindCell()
                        {
                                iAx = 1, 0, 0; 
                                iAy = 0, 1, 0; 
                                iAz = 0, 0, 1;
                        }
                
                virtual ~FindCell(){}
                
                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 ) = 0;
                virtual point  localUVW   ( const point&p ) = 0;
        };

        
        class HexaHedralCellNewton : public FindCell
        {
                RefPtr<UniGridMapper> UniMap;

                point uvw;
                int point_type;
                double block_epsilon;

        public:
                HexaHedralCellNewton():FindCell()
                {}
 
                virtual ~HexaHedralCellNewton(){}

                void  setUniMap  ( RefPtr<UniGridMapper> unimap ) { UniMap = unimap; }
                int   pointInCell( const point&P );
                point localUVW   ( const point&p );

                void setEpsilon(const double e) { block_epsilon = e; }
        }; 


        bool rayTriangleIntersect(  const point orig, const tvector dir, 
                                    const point vert0, const point vert1, const point vert2, 
                                    double&t, double&u, double&v,
                                    point&intersection );
                

        bool rayQuaderIntersect(  const point orig, const tvector dir, 
                                  Vertexlist_t & QV,
                                  double&t, double&u, double&v,
                                  point&intersection, bivector&normal, unsigned&face ); 

        FixedArray<FixedArray<int, 3>,12> Facelist; 

        // to be removed later
        double block_epsilon; 

                // parameters for unigrid building
        double  res_scale, prec_scale;

public: 
        int     MaxListSize;

        LocalFromWorldPoint( Slice&SourceSliceP, const RefPtr<Grid>&SourceGridP,
                             const string&Gridname,
                             const RefPtr<Field>&SourceCoordsP,
                             double     res_scale = 1.0, double prec_scale = 0.1);

        ~LocalFromWorldPoint();

        bool get( const point position, pair<point, string>&data ); 


        /* Set an epsilon value for widening the block boundary test. This is a hack solution and should be done internally via
         * neighbourhood information avoiding espsilons. Some tests showed that a value of 0.05 seems to do pretty ok.
         */
        void setBlockBoundaryEpsilon(const double block_epsilonP) { block_epsilon = block_epsilonP; }


        int numberOfBordersOfCell( const pair<point, string> localdata, int&hint1, int&hint2 ); 


        bool getBorderIntersectionNormal( const point worldposition, 
                                          const pair<point, string> localpoint, 
                                          const tvector direction,  
                                          point&intersection, bivector&normal ); 
};

} // namespace

#endif