GridEvaluator_impl.hpp

#ifndef __FIBER_BASEOPERATIONS_GRIDEVALUATOR_IMPL_HPP
#define __FIBER_BASEOPERATIONS_GRIDEVALUATOR_IMPL_HPP
 
#include "GridEvaluator.hpp"
#include <field/DirectProductArray.hpp>
#include <field/ArrayRef.hpp>
 
#include <grid/types/FEM.hpp>
 
#include <vector/Interpolate.hpp>
#include <vector/LinearIpol.hpp>
 
#include "LocalFromWorldPoint.hpp"
 
#include <aerie/KDTree.hpp>
 
#include <field/FragmentIDContainer.hpp>
 
namespace Fiber
{
 
struct  FragmentSearchBase
{
        GridEvaluator::Context&myContext;
        RefPtr<Field>   SourceCoordinates;
 
        RefPtr<LocalFromWorldPoint>     myLocalFromWorldPoint;
 
 
        RefPtr<FragmentID>      FoundFragID;
 
        FragmentSearchBase(GridEvaluator::Context&theContext)
        : myContext(theContext)
        {}
 
        LocalFromWorldPoint&getLocalFromWorldPoint()
        {
                if (!myLocalFromWorldPoint)
                {
                        myLocalFromWorldPoint = new LocalFromWorldPoint(myContext.mySlice, myContext.mySourceGrid, myContext.Gridname
//                                                                      , SourceCoordinates
                                                                        , myContext.Res_scale, myContext.Prec_scale);
                }
                return *myLocalFromWorldPoint;
        }
 
 
};
 
template <class PointType>
        struct  CurvilinearSearch
{
 
template <class FloatIndexType>
static  bool    getPoint(FragmentSearchBase&FSB ,
                         const PointType&P,
                         FloatIndexType&InterpolationWeights,
                         RefPtr<FragmentID>&FragID)
        {
                return false;
        }
};
 
template <>
struct  CurvilinearSearch<Eagle::PhysicalSpace::point>
{
        typedef FixedArray<double, 3> FloatIndex_t;
 
static  bool    getPoint(FragmentSearchBase&FSB,
                         const Eagle::PhysicalSpace::point&P,
                         FloatIndex_t&InterpolationWeights,
                         RefPtr<FragmentID>&FragID)
        {
        LocalFromWorldPoint&findLocal = FSB.getLocalFromWorldPoint(); 
 
//      pair<Eagle::PhysicalSpace::point, RefPtr<FragmentID> > MultiblockFloatingIndex; 
        LocalPoint MultiblockFloatingIndex;
        bool    FoundAndEvaluatedThePoint = findLocal.get(P, MultiblockFloatingIndex); 
 
                InterpolationWeights = MultiblockFloatingIndex.local_point; 
                FragID               = MultiblockFloatingIndex.frag_id; 
                return FoundAndEvaluatedThePoint;
        }
};
 
 
// some structure for KD search, maybe
struct VertexWeights
{ 
        std::vector<index_t> index; 
        std::vector<double>  weight; 
};
 
template <class PointType, class FloatIndexType>
struct KDTreeCloudSearch
{
 
 
        static  bool    getPoint( RefPtr<MemArray<1, PointType> >&VertexSet,
                         const PointType&P,
                         FloatIndexType&InterpolationWeights,
                         RefPtr<FragmentID>&FragID)
        {
                puts("KDTreeCloudSearch: wrong weight type requested!!\n"); 
                return false;
        }
};
 
template <class PointType>
struct KDTreeCloudSearch<PointType, VertexWeights>
{
enum { N = PointType::SIZE };
 
        static  bool    getPoint( RefPtr<MemArray<1, PointType> >&VertexSet,
                         const PointType&P,
                         VertexWeights&InterpolationWeights,
                         RefPtr<FragmentID>&FragID)
        {
                // Do the search here ...
 
                RefPtr<Eagle::KDTree<N, index_t> > myTree ; 
                RefPtr<Eagle::KDInterface <N, index_t> > TreeInt = VertexSet->findInterface( typeid (Eagle::KDTree<2, index_t>) ) ; 
                        // check for a cached tree on the vertexset. fill and cache tree if not found 
                        if( !TreeInt )
                        {
#ifdef VERBOSE
                                puts("GridEvaluator::creating new tree"); 
#endif 
                                myTree = new Eagle::KDTree<N, index_t >() ; 
                                TreeInt = new Eagle::KDInterface <N, index_t >( myTree ) ; 
 
                        MultiArray<1, PointType> Vertices = *VertexSet; 
                        MultiIndex<1> max = Vertices.Size(); 
                        MultiIndex<1> mi(0);
 
                                do
                                {
                                        myTree->insert( Vertices[mi], mi[0] );
                                }while ( mi.inc(max) );
 
                                VertexSet->addInterface( TreeInt ) ; 
                        } 
                        else // tree found -> assign it to local tree pointer
                        {
#ifdef VERBOSE
                                puts("GridEvaluator::found cached tree"); 
#endif 
                                myTree = TreeInt->Tree ; 
                        } 
 
                // data container for result query poitns 
                std::multimap<double, index_t>verts; 
 
                // somehow a radius for the query must be controlled, or also be adjustable paramter.. 
                double r = 1.0;
 
                        // check for the certain poitn P 
                        myTree->getEuclideanRange(P, r, verts ); 
 
                        // no candidates found in range
                        if( verts.size() < 0 )
                                return GridEvaluator::Nope; 
 
                        // copy weights into VertexWeightstructure 
                        // maybe better implement container filling directly in KDTree 
                       // (but then it would be unsorted ... )
                index_t i = 0;
                        InterpolationWeights.weight.resize( verts.size() );
                        InterpolationWeights.index.resize( verts.size() );
 
                        for ( std::multimap<double, index_t>::iterator it = verts.begin(); it != verts.end(); it++ )
                        {
                                InterpolationWeights.weight[i] = (*it).first; 
                                InterpolationWeights.index[i]  = (*it).second; 
                                i++;
                        }
 
                return true;
        }
};
 
 
 
// mr: floatindex_t should become a template parameter, such it can be exchanged be a more general type?
template <class PointType>
struct  FragmentSearch : FragmentSearchBase
{
        enum { N = PointType::SIZE };
 
        typedef FixedArray<double, N> FloatIndex_t;
 
 
        FragmentSearch(GridEvaluator::Context&theContext)
        : FragmentSearchBase(theContext)
        {}
 
 
        GridEvaluator::FragmentSearchResult Contains(const PointType&P,
                                                     const MemCore::StrongPtr<Fiber::CreativeArrayBase>&DC,
                                                     FloatIndex_t&FloatIndex)
        {
                // try to get data, if loaded.
        RefPtr<MemBase> AttemptedFragmentCoordinates = DC->get();
 
                // TODO: Check for Min/Max Interface stored at DC
 
        RefPtr<MemBase> FragmentCoordinates = DC->create(); // possibly slow call! 
                if (!FragmentCoordinates)
                {
                        assert( 0 && "bah!");
                        return GridEvaluator::Nope; 
                }
 
                if (RefPtr<DirectProductMemArray<PointType, LinearCoordinateAxisArray> >
                        UniformCoordinates = FragmentCoordinates)
                {
                        FloatIndex = UniformCoordinates->getIndex( P ); 
/*
                        for(int i=0; i<N; i++)
                        {
                                if (FloatIndex[i]<0.0 || FloatIndex[i]>1.0)
                                        return GridEvaluator::Nope;
                        } 
*/
                        return GridEvaluator::YesWithLocationAlreadyFound;
                } 
 
                if (RefPtr<DirectProductMemArray<PointType, ExplicitCoordinateAxisArray> >
                        UniformCoordinates = FragmentCoordinates)
                {
                        FloatIndex = UniformCoordinates->getIndex( P ); 
 
                        for(int i=0; i<N; i++)
                        {
                                if (FloatIndex[i]<0.0 || FloatIndex[i]>1.0)
                                        return GridEvaluator::Nope;
                        }
 
                        return GridEvaluator::YesWithLocationAlreadyFound;
                } 
 
                if(  RefPtr<MemArray<1, PointType> > VertexSet = FragmentCoordinates ) // 1D explicit vertices
                {
//                      if( (*myContext.mySourceGrid)(FEM::cellID()) 
//                          && (*(*myContext.mySourceGrid)(FEM::cellID() ))( myContext.mySourceGrid->findVertices() ) // have a relative representation cells as vertices
//                          && ((*(*myContext.mySourceGrid)(FEM::cellID() ))( myContext.mySourceGrid->findVertices() ))->Positions() // have cell positions field
//                          && ((*(*myContext.mySourceGrid)(FEM::cellID() ))( myContext.mySourceGrid->findVertices() ))->Positions()->getData( FoundFragID ) // have data
//                             // support for fragmented data? same frags in cells as in verts
//                              )
//                      {
//                              // build/check for search structures: vertices as cells -> converter? cell neighborhood? -> converter?
//                              // max cell size, cell midpoints? KD query?
//                              // maybe remove mixed type meshes and assume one element type per part (which is correct for abaqus) -> converter?
//                              // weight would be similar to the point cloud case with indices and weights... but maybe not vlens but fixed sizes
                                
//                      RefPtr<MemArray<1, CellIndices_t> > cell_set = ((*(*myContext.mySourceGrid)(FEM::cellID() ))( myContext.mySourceGrid->findVertices() ))->Positions()->getData( FoundFragID ); 
 
// //                   RefPtr<Eagle::KDTree<PointType::SIZE, index_t> > tree = FEM::getCellMidPointTree<PointType>( cell_set, VertexSet ); 
// //                   double max_cell_size = 1.0; //FEM::maxCellSize( cell_set ); 
 
// //                   std::list<index_t> candidates; 
// //                           tree->nearest_range(P, max_cell_size, candidates ); 
 
// //                   std::list<index_t>::const_iterator it; 
// //                           for( it = candidates.begin(); it != candidates.end(); it++ ) 
// //                           {
//                                      // do some check to find local coords 
// //                           }
//                      } 
//                      // else ifs: check other special 1D explicit point3D types 
//                      else // point cloud
//                      {
                                puts("GridEvaluator: Evaluating set of Vertices not yet fully implemented!!!\n"); 
                                fflush(stdout); 
 
                                // mr: so what to come to the VertexWeighting type here?
                                if (KDTreeCloudSearch<PointType, FloatIndex_t >::getPoint( VertexSet, P, FloatIndex, FoundFragID) )
                                {
//                              return GridEvaluator::YesWithLocationAlreadyFound; 
                                        return  GridEvaluator::Nope;
                                } 
 
 
                                Verbose(0) << "###### GridEvaluator: Evaluating set of Vertices not yet implemented!!! ####";
                                return GridEvaluator::Nope; 
//                      }
 
                }
 
                /*TODO: Marcel, this is the place where you use the KD-tree lookup to find
                  all points in the VertexSet which are close to the provided point "P".
                  Currently the result is expected to be returned in a FloatIndex_t . 
                  This will NOT be appropriate for an evaluation of a point set and
                  we will need to extend this API for a more general evaluation weight,
                  but need careful to not impact performance. In the most general case
                  the interpolation weight will be an index_t (integer) with a distance
                  weight (double), thus replacing FloatIndex_t with
                  struct { vector<index_t> index; vector<double> weight; };
                  but this semantic structure is not very efficient due to the
                  dynamic allocation of vector's. Need to think about a more efficient
                  scheme.
                 */ 
/*
                // mr: Couldnt this also be a line instead of a pointcloud?
                if (RefPtr<MemArray<1, PointType> > VertexSet = FragmentCoordinates)
                {
                        puts("GridEvaluator: Evaluating set of Vertices not yet fullly implemented!!!\n"); 
                        fflush(stdout); 
 
                        // mr: to be adjusted to apropriate weight type
                        if (KDTreeCloudSearch<PointType, FloatIndex_t >::getPoint( VertexSet, P, FloatIndex, FoundFragID) )
                        {
//                              return GridEvaluator::YesWithLocationAlreadyFound; 
                                return  GridEvaluator::Nope;
                        } 
 
 
                        assert( 0 && "GridEvaluator: Evaluating set of Vertices not yet implemented!!!\n"); 
                        return GridEvaluator::Nope;
                } 
*/
                if (RefPtr<MemArray<3, PointType> > CurviLinear3DCoords = FragmentCoordinates)
                {
 
                        if (CurvilinearSearch<PointType>::getPoint( *this, P, FloatIndex, FoundFragID) )
                        {
                                puts("yayy, curvigrid found!"); 
                                return GridEvaluator::YesWithLocationAlreadyFound;
                        } 
 
                        puts("GridEvaluator: NO curvilinear 3D location found!\n"); 
//                      fflush(stdout); 
//                      assert( 0 && "GridEvaluator: Evaluating curvilinear 3D not yet implemented!!!\n");
                        return GridEvaluator::Nope;
                }
 
 
                // curvilinear? do coarse bounding box check.
 
                assert( 0 && "GridEvaluator: Provided Coordinate type not supported!!!\n"); 
                return GridEvaluator::Nope;
        }
};
 
 
 
 
 
template <class PositionsFieldType>
void    GridEvaluator::findForSpecificDestinationPositions(const PositionsFieldType&DestinationPositions, 
                                                           const Field&SourcePositions,
                                                           GridEvaluator::Context&theEvaluationContext)
{
        typedef typename PositionsFieldType::value_type point_t; 
 
        // the dimensionality of the destination manifold, for instance 2 for a plane
        enum { DestManifoldDims      = PositionsFieldType::Dims }; 
 
        // the dimensionality of the manifold where this plane is embedded, for instance 
        // 3 for a plane in a 3D space.
        enum { EmbeddingManifoldDims = point_t::SIZE }; 
 
typedef FixedArray<double, EmbeddingManifoldDims> EmbeddingFloatIndex_t; 
 
typedef MemArray  <DestManifoldDims, FragmentIDContainer> FragmentMemArray_t;
typedef MultiArray<DestManifoldDims, FragmentIDContainer> FragmentMultiArray_t;
 
MultiIndex<DestManifoldDims> NumberOfEvalPositions = DestinationPositions.Size(); 
 
        if (NumberOfEvalPositions.size()<1)
                return;
 
MemCore::Ref<FragmentMemArray_t> FragmentsPerPosition( NumberOfEvalPositions );
FragmentMultiArray_t&FragmentArray = *FragmentsPerPosition; 
 
ArrayRef<EmbeddingFloatIndex_t, DestManifoldDims> InterpolationWeightArray( NumberOfEvalPositions );
ArrayRef<int, DestManifoldDims> FragmentSearchStatus( NumberOfEvalPositions );
 
 
        /*
          Iterate over all fragments of the source field
          to collect all possible fragment ID's for a given 
          destination point.
         */
        struct  SourceFragmentIterator : Field::Iterator
        {
                FragmentSearch<point_t> FragmentSearcher;
 
                const point_t        &DestPoint;
                FragmentIDContainer  &SourceFragments;
                int                  &SearchStatus;
                FixedArray<double, EmbeddingManifoldDims>&DestIndexInSource;
 
                SourceFragmentIterator( GridEvaluator::Context&theEvaluationContext,
                                        const point_t        &theDestPoint,
                                        FragmentIDContainer  &theSourceFragments,
                                        int                  &theSearchStatus,
                                        FixedArray<double, EmbeddingManifoldDims>&theDestIndexInSource)
                : FragmentSearcher( theEvaluationContext )
                , DestPoint( theDestPoint )
                , SourceFragments( theSourceFragments )
                , SearchStatus ( theSearchStatus )
                , DestIndexInSource( theDestIndexInSource )
                {
                        SearchStatus = Nope;
                }
 
                bool apply(const RefPtr<FragmentID>&f, const MemCore::StrongPtr<Fiber::CreativeArrayBase>&DC)
                {
                        if (!DC) return true; 
 
                FragmentSearchResult GotIt = FragmentSearcher.Contains( DestPoint, DC, DestIndexInSource ); 
 
                        SearchStatus = GotIt;
 
                        switch( GotIt )
                        {
                        case Possibly:
                        case Certainly:
                        case YesWithLocationAlreadyFound:
                                
                                if (FragmentSearcher.FoundFragID)
                                {
                                        SourceFragments.insert(FragmentSearcher.FoundFragID);
                                        return false; // stop iteration over fragments for the destination point
                                }
 
                                SourceFragments.insert( f );
                                break;
 
                        default: ;
                        }
 
                        return true;
                }
        };
 
// Iteration over all *destination* points 
// @TODO should be done in parallel! 
// @TODO need to optimize to move out of this loop as much as posisble
MultiIndex<DestManifoldDims> Index;
        do
        {
        int                                      &theSearchStatus      = FragmentSearchStatus[ Index ]; 
        FixedArray<double, EmbeddingManifoldDims> theDestIndexInSource;
 
        // we want to know the fragment ID's at this point 
        SourceFragmentIterator GetFragmentsForPoint( theEvaluationContext,
                                                     DestinationPositions[ Index ],
                                                     FragmentArray[ Index ],
 
                                                     theSearchStatus,
                                                     theDestIndexInSource ); 
 
                GetFragmentsForPoint.FragmentSearcher.SourceCoordinates = SourcePositions.self();
 
                SourcePositions.iterate( GetFragmentsForPoint );
 
                InterpolationWeightArray[ Index ] = theDestIndexInSource; 
 
//              printf("handled index [%d,%d] stat %d\n", Index[0], Index[1], theSearchStatus);
        }
        while( Index.inc( NumberOfEvalPositions ) ); 
 
 
        FragmentIDs   = FragmentsPerPosition; 
        SearchStatus  = FragmentSearchStatus;
        PointLocation = InterpolationWeightArray;
}
 
 
template <int SourceDims, class ValueType>
bool    GridEvaluator::evalFloatIndices(std::vector<ValueType>&ResultData,
                                        const CreativeIterator< FixedArray<double, SourceDims> >&InterpolationWeightsPerDestPoint,
                                        const Field&SourceField) const
{
typedef FixedArray<double, SourceDims> FloatIndex_t;
 
TypedIterator<FragmentIDContainer> FragmentsPerDestPoint  = FragmentIDs; 
 
        assert(FragmentsPerDestPoint);
 
        if (TypedIterator<int> SearchResult = SearchStatus)
        {
        const CreativeIterator<int>&SR = *SearchResult;
        const Iterator<int>&It = SR; 
 
        const CreativeIterator<FragmentIDContainer>&FragmentPerDestPoint = *FragmentsPerDestPoint;
 
/*
  A different strategy would be to first collect all destination points
  that belong to the same source fragment and then evaluate all those
  points from the same source fragment. This might be more performant
  because the fragment needs to be retrieved only once instead of for
  each point. It would also be more memory-efficient if only
  one fragment can be kept in memory.
 */ 
 
                for(index_t i=0; i<It.count(); i++)
                {
                        switch( FragmentSearchResult(It[ i ] ) )
                        {
                        case YesWithLocationAlreadyFound:
                                {
                                const FloatIndex_t       &PointIndex     = InterpolationWeightsPerDestPoint[ i ]; 
                                const FragmentIDContainer&SourceFragment = FragmentPerDestPoint[i]; 
 
                                        // TODO: reverse loop such that this can be outside, then do in parallel
                                        if (SourceFragment.nFragments()==1)
                                        {
                                        RefPtr<FragmentID> fID = SourceFragment.getFirst();
                                                // note: fID can be a NullPtr, this is valid, and indicates
                                                // an unfragmented source field. 
                                        RefPtr<MemArray<SourceDims, ValueType> > 
                                                SourceMemArray = SourceField.getData( fID ); 
                                                if (!SourceMemArray) return false;
 
                                        MultiArray<SourceDims, ValueType>  SourceData = *SourceMemArray; 
 
//                                      Interpolate<SourceDims, ValueType, LinearIpol<ValueType>, ValueType> 
                                        Interpolate<SourceDims, ValueType, LinearIpol<ValueType> >
                                                MyInterpolator(SourceData, PointIndex);
 
                                        const ValueType&InterpolatedValue = MyInterpolator.eval(); 
                                                ResultData[ i ] = InterpolatedValue; 
 
                                        } 
                                        assert( SourceFragment.nFragments() == 1); 
                                        break;
                                }
 
                        default:
                                //              assert(0 && "Uncertain fragment not handled yet");
                                break;
                                ;
                        }
                } 
 
                return true;
        }
 
        assert( 0 && "No point fragment location per point available, not good.");
 
        {
        }
 
        // type: FixedArray<double, int CoordinateDims> 
//      RefPtr<MemBase> PointLocation; 
 
        return false;
}
 
 
 
template <class ValueType>
bool    GridEvaluator::evalType(std::vector<ValueType>&ResultData, const Field&SourceField) const
{
        if (!PointLocation)
                return false; 
 
        if ( TypedIterator< FixedArray<double, 3> > EmbeddingManifoldIs3D = PointLocation) 
        {
                return evalFloatIndices(ResultData, *EmbeddingManifoldIs3D, SourceField);
        } 
        if ( TypedIterator< FixedArray<double, 2> > EmbeddingManifoldIs2D = PointLocation) 
        {
                return evalFloatIndices(ResultData, *EmbeddingManifoldIs2D, SourceField);
        }
 
        return false;
}
 
 
 
} // namespace Fiber
 
 
#endif // __FIBER_BASEOPERATIONS_GRIDEVALUATOR_IMPL_HPP