hit_test_factory.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/geometry/target_set.hpp>
#include <kernel/geometry/index_set.hpp>
#include <kernel/geometry/mesh_part.hpp>
#include <kernel/geometry/atlas/chart.hpp>
#include <kernel/util/tiny_algebra.hpp>
 
namespace FEAT
{
  namespace Geometry
  {
    namespace Intern
    {
      template<typename VertexSet_, typename IndexTuple_>
      inline Tiny::Vector<typename VertexSet_::CoordType, VertexSet_::num_coords>
      get_midpoint(const VertexSet_& vertex_set, const IndexTuple_& index_tuple)
      {
        using DataType = typename VertexSet_::CoordType;
        constexpr int num_vertex = IndexTuple_::num_indices;
 
        Tiny::Vector<DataType, VertexSet_::num_coords> mid_point(DataType(0));
        for(int i(0); i < num_vertex; ++i)
        {
          for(int j(0); j < VertexSet_::num_coords; ++j)
          {
            mid_point[j] += vertex_set[index_tuple[i]][j];
          }
        }
        return mid_point * (DataType(1) / DataType(num_vertex));
      }
 
      template<typename VertexSet_>
      static Tiny::Vector<typename VertexSet_::CoordType, VertexSet_::num_coords>
      get_midpoint(const VertexSet_& vertex_set, const Index vertex)
      {
        using DataType = typename VertexSet_::CoordType;
 
        Tiny::Vector<DataType, VertexSet_::num_coords> mid_point;
        for(int j(0); j < VertexSet_::num_coords; ++j)
        {
          mid_point[j] = vertex_set[vertex][j];
        }
        return mid_point;
      }
 
      template<typename HitFunc_, typename Mesh_, int dim_ = Mesh_::shape_dim>
      inline void hittest_compute_target_data(
        std::vector<std::vector<Index>>& target_data,
        const Mesh_& mesh,
        const HitFunc_& hit_test)
      {
        if constexpr (dim_ == 0)
        {
          const Index num_cells(mesh.get_num_entities(0));
 
          target_data[0].reserve(num_cells);
          for(Index i(0); i < num_cells; ++i)
          {
            if (hit_test(get_midpoint(mesh.get_vertex_set(), i)))
            {
              target_data[0].push_back(i);
            }
          }
        }
        else
        {
          // Recurse
          hittest_compute_target_data<HitFunc_, Mesh_, dim_ - 1>(target_data, mesh, hit_test);
 
          const Index num_cells(mesh.get_num_entities(dim_));
          const auto& index_set(mesh.template get_index_set<dim_,0>());
 
          target_data[dim_].reserve(num_cells);
          for(Index i(0); i < num_cells; ++i)
          {
            if (hit_test(get_midpoint(mesh.get_vertex_set(), index_set[i])))
            {
              target_data[dim_].push_back(i);
            }
          }
        }
      }
 
      template<typename HitFunc_, typename Mesh_, int dim_ = Mesh_::shape_dim>
      inline void hittest_compute_filtered_target_data(
        std::vector<std::vector<Index>>& target_data,
        const Mesh_& mesh,
        const MeshPart<Mesh_>& part,
        const HitFunc_& hit_test)
      {
        if constexpr (dim_ == 0)
        {
          const Index num_cells(part.get_num_entities(0));
          const TargetSet& vertex_target_set = part.template get_target_set<0>();
 
          target_data[0].reserve(num_cells);
          for(Index i(0); i < num_cells; ++i)
          {
            const Index vertex = vertex_target_set[i];
            if (hit_test(get_midpoint(mesh.get_vertex_set(), vertex)))
            {
              target_data[0].push_back(vertex);
            }
          }
        }
        else
        {
          // Recurse
          hittest_compute_filtered_target_data<HitFunc_, Mesh_, dim_ - 1>(target_data, mesh, part, hit_test);
 
          const Index num_cells(part.get_num_entities(dim_));
          const TargetSet& filter_target_set = part.template get_target_set<dim_>();
          const auto& index_set(mesh.template get_index_set<dim_,0>());
 
          target_data[dim_].reserve(num_cells);
          for(Index i(0); i < num_cells; ++i)
          {
            const Index entity = filter_target_set[i];
            if (hit_test(get_midpoint(mesh.get_vertex_set(), index_set[entity])))
            {
              target_data[dim_].push_back(entity);
            }
          }
        }
      }
 
      template<typename Shape_, int dim_ = Shape_::dimension>
      inline void write_to_target_set(TargetSetHolder<Shape_>& tsh, std::vector<std::vector<Index>>& source)
      {
        if constexpr (dim_ < 0)
        {
          return;
        }
        else
        {
          // Recurse
          write_to_target_set<Shape_, dim_ - 1>(tsh, source);
 
          const auto num_cells(Index(source[dim_].size()));
          TargetSet& trg = tsh.template get_target_set<dim_>();
          for(Index i(0); i < num_cells; ++i)
          {
            trg[i] = source[dim_][i];
          }
 
        }
      }
    } // namespace Intern
 
    template<
      typename HitFunc_,
      typename Mesh_>
    class HitTestFactory :
      public Factory< MeshPart<Mesh_> >
    {
    public:
      typedef typename Mesh_::ShapeType ShapeType;
      typedef MeshPart<Mesh_> MeshType;
      typedef typename MeshType::TargetSetHolderType TargetSetHolderType;
 
    protected:
      const HitFunc_& _hit_func;
      const Mesh_& _mesh;
      std::vector<std::vector<Index>> _target_data;
 
    public:
      explicit HitTestFactory(const HitFunc_& hit_func, const Mesh_& mesh) :
        _hit_func(hit_func),
        _mesh(mesh),
        _target_data(std::size_t(_mesh.shape_dim+1))
      {
        Intern::hittest_compute_target_data<HitFunc_, Mesh_>(_target_data, _mesh, _hit_func);
      }
 
      explicit HitTestFactory(const HitFunc_& hit_func, const Mesh_& mesh, const MeshType& filter) :
        _hit_func(hit_func),
        _mesh(mesh),
        _target_data(std::size_t(_mesh.shape_dim+1))
      {
        Intern::hittest_compute_filtered_target_data<HitFunc_, Mesh_>(_target_data, _mesh, filter, _hit_func);
      }
 
      virtual Index get_num_entities(int dim) override
      {
        return Index(_target_data.at(std::size_t(dim)).size());
      }
 
      virtual void fill_target_sets(TargetSetHolderType& target_set_holder) override
      {
        Intern::write_to_target_set(target_set_holder, _target_data);
      }
 
      virtual void fill_attribute_sets(typename MeshType::AttributeSetContainer&) override
      {
        // do nothing as the object has no attribute sets
      }
 
      virtual void fill_index_sets(std::unique_ptr<typename MeshType::IndexSetHolderType>&) override
      {
        // do nothing as the object has no index sets
      }
 
    }; // class HitTestFactory
 
    template<typename Mesh_, typename HitFunc_>
    MeshPart<Mesh_> make_meshpart_by_hit_test(const Mesh_& mesh, const HitFunc_& hit_func)
    {
      HitTestFactory<HitFunc_, Mesh_> factory(mesh, hit_func);
      return factory.make();
    }
 
    template<typename Mesh_, typename HitFunc_>
    MeshPart<Mesh_> make_meshpart_by_filtered_hit_test(const Mesh_& mesh, const MeshPart<Mesh_>& filter, const HitFunc_& hit_func)
    {
      HitTestFactory<HitFunc_, Mesh_> factory(mesh, filter, hit_func);
      return factory.make();
    }
 
    template<typename Mesh_, typename HitFunc_>
    std::unique_ptr<MeshPart<Mesh_>> make_unique_meshpart_by_hit_test(const Mesh_& mesh, const HitFunc_& hit_func)
    {
      HitTestFactory<HitFunc_, Mesh_> factory(mesh, hit_func);
      return factory.make_unique();
    }
 
    template<typename Mesh_, typename HitFunc_>
    std::unique_ptr<MeshPart<Mesh_>> make_unique_meshpart_by_filtered_hit_test(const Mesh_& mesh, const MeshPart<Mesh_>& filter, const HitFunc_& hit_func)
    {
      HitTestFactory<HitFunc_, Mesh_> factory(mesh, filter, hit_func);
      return factory.make_unique();
    }
 
    template<typename DataType_, int dim_>
    class SphereHitTestFunction
    {
    public:
      typedef Tiny::Vector<DataType_, dim_> PointType;
 
    private:
      PointType _midpoint;
      DataType_ _radius;
 
    public:
      explicit SphereHitTestFunction(DataType_ radius) :
        _midpoint(DataType_(0)),
        _radius(radius)
      {
        XASSERTM(radius > DataType_(0), "sphere radius must be > 0");
      }
 
      SphereHitTestFunction(PointType midpoint, DataType_ radius) :
        _midpoint(midpoint),
        _radius(radius)
      {
        XASSERTM(radius > DataType_(0), "sphere radius must be > 0");
      }
 
      bool operator()(PointType point) const
      {
        return (point - _midpoint).norm_euclid() <= _radius;
      }
    }; // class SphereHitTestFunction<...>
 
 
    template<typename Mesh_, typename Chart_ = Atlas::ChartBase<Mesh_>>
    class ChartHitTestFactory :
      public Factory< MeshPart<Mesh_> >
    {
    public:
      typedef typename Mesh_::ShapeType ShapeType;
      typedef MeshPart<Mesh_> MeshType;
      typedef typename MeshType::TargetSetHolderType TargetSetHolderType;
 
      class ChartHitFunction
      {
        typedef typename Mesh_::CoordType CoordType;
        typedef typename Mesh_::VertexType WorldPointType;
 
      protected:
        const Geometry::Atlas::ChartBase<Mesh_>& _chart;
        const CoordType _inv;
        const CoordType _tol;
 
      public:
        explicit ChartHitFunction(const Geometry::Atlas::ChartBase<Mesh_>& chart, bool invert, CoordType tol = CoordType(0)) :
          _chart(chart),
          _inv(CoordType(invert ? -1 : 1)),
          _tol(tol)
        {
        }
 
        bool operator()(const WorldPointType& point) const
        {
          return _inv * _chart.signed_dist(point) <= _tol;
        }
      };
 
    protected:
      const Mesh_& _mesh;
      ChartHitFunction _hit_func;
      std::vector<std::vector<Index>> _target_data;
 
    public:
      explicit ChartHitTestFactory(const Mesh_& mesh, const Chart_& chart, bool invert = false, typename Mesh_::CoordType tol = typename Mesh_::CoordType(0)) :
        _mesh(mesh),
        _hit_func(chart, invert, tol),
        _target_data(std::size_t(_mesh.shape_dim+1))
      {
        Intern::hittest_compute_target_data<ChartHitFunction, Mesh_>(_target_data, _mesh, _hit_func);
      }
 
      explicit ChartHitTestFactory(const Mesh_& mesh, const MeshType& filter, const Chart_& chart, bool invert = false, typename Mesh_::CoordType tol = typename Mesh_::CoordType(0)) :
        _mesh(mesh),
        _hit_func(chart, invert, tol),
        _target_data(std::size_t(_mesh.shape_dim+1))
      {
        Intern::hittest_compute_filtered_target_data<ChartHitFunction, Mesh_>(_target_data, _mesh, filter, _hit_func);
      }
 
      virtual Index get_num_entities(int dim) override
      {
        return Index(_target_data.at(std::size_t(dim)).size());
      }
 
      virtual void fill_target_sets(TargetSetHolderType& target_set_holder) override
      {
        Intern::write_to_target_set<ShapeType>(target_set_holder, _target_data);
      }
 
      virtual void fill_attribute_sets(typename MeshType::AttributeSetContainer&) override
      {
        // do nothing as the object has no attribute sets
      }
 
      virtual void fill_index_sets(std::unique_ptr<typename MeshType::IndexSetHolderType>&) override
      {
        // do nothing as the object has no index sets
      }
    }; // class ChartHitTestFactory
 
    template<typename Mesh_, typename Chart_>
    MeshPart<Mesh_> make_meshpart_by_chart_hit_test(const Mesh_& mesh, const Chart_& chart, bool invert = false, typename Mesh_::CoordType tol = typename Mesh_::CoordType(0))
    {
      ChartHitTestFactory<Mesh_, Chart_> factory(mesh, chart, invert, tol);
      return factory.make();
    }
 
    template<typename Mesh_, typename Chart_>
    MeshPart<Mesh_> make_meshpart_by_filtered_chart_hit_test(const Mesh_& mesh, const MeshPart<Mesh_>& filter, const Chart_& chart, bool invert = false, typename Mesh_::CoordType tol = typename Mesh_::CoordType(0))
    {
      ChartHitTestFactory<Mesh_, Chart_> factory(mesh, filter, chart, invert, tol);
      return factory.make();
    }
 
    template<typename Mesh_, typename Chart_>
    std::unique_ptr<MeshPart<Mesh_>> make_unique_meshpart_by_chart_hit_test(const Mesh_& mesh, const Chart_& chart, bool invert = false, typename Mesh_::CoordType tol = typename Mesh_::CoordType(0))
    {
      ChartHitTestFactory<Mesh_, Chart_> factory(mesh, chart, invert, tol);
      return factory.make_unique();
    }
 
    template<typename Mesh_, typename Chart_>
    std::unique_ptr<MeshPart<Mesh_>> make_unique_meshpart_by_filtered_chart_hit_test(const Mesh_& mesh, const MeshPart<Mesh_>& filter, const Chart_& chart, bool invert = false, typename Mesh_::CoordType tol = typename Mesh_::CoordType(0))
    {
      ChartHitTestFactory<Mesh_, Chart_> factory(mesh, filter, chart, invert, tol);
      return factory.make_unique();
    }
 
    template<typename Mesh_, typename Lambda_>
    class LambdaHitTestFactory :
      public Factory< MeshPart<Mesh_> >
    {
    public:
      typedef typename Mesh_::ShapeType ShapeType;
      typedef MeshPart<Mesh_> MeshType;
      typedef typename MeshType::TargetSetHolderType TargetSetHolderType;
 
    protected:
      class HitTestLambda
      {
      public:
        Lambda_ lambda;
 
        explicit HitTestLambda(Lambda_&& lb) :
          lambda(std::forward<Lambda_>(lb))
        {
        }
 
        template<typename T_, int s_>
        bool operator()(const Tiny::Vector<T_, 1, s_>& p) const
        {
          return lambda(p[0]);
        }
 
        template<typename T_, int s_>
        bool operator()(const Tiny::Vector<T_, 2, s_>& p) const
        {
          return lambda(p[0], p[1]);
        }
 
        template<typename T_, int s_>
        bool operator()(const Tiny::Vector<T_, 3, s_>& p) const
        {
          return lambda(p[0], p[1], p[3]);
        }
      } _hit_func;
      const Mesh_& _mesh;
      std::vector<std::vector<Index>> _target_data;
 
    public:
      explicit LambdaHitTestFactory(const Mesh_& mesh,  Lambda_&& lambda) :
        _hit_func(std::forward<Lambda_>(lambda)),
        _mesh(mesh),
        _target_data(std::size_t(_mesh.shape_dim+1))
      {
        Intern::hittest_compute_target_data<HitTestLambda, Mesh_>(_target_data, _mesh, _hit_func);
      }
 
      explicit LambdaHitTestFactory(const Mesh_& mesh, const MeshType& filter, Lambda_&& lambda) :
        _hit_func(std::forward<Lambda_>(lambda)),
        _mesh(mesh),
        _target_data(std::size_t(_mesh.shape_dim+1))
      {
        Intern::hittest_compute_filtered_target_data<HitTestLambda, Mesh_>(_target_data, _mesh, filter, _hit_func);
      }
 
      virtual Index get_num_entities(int dim) override
      {
        return Index(_target_data.at(std::size_t(dim)).size());
      }
 
      virtual void fill_target_sets(TargetSetHolderType& target_set_holder) override
      {
        Intern::write_to_target_set<ShapeType>(target_set_holder, _target_data);
      }
 
      virtual void fill_attribute_sets(typename MeshType::AttributeSetContainer&) override
      {
        // do nothing as the object has no attribute sets
      }
 
      virtual void fill_index_sets(std::unique_ptr<typename MeshType::IndexSetHolderType>&) override
      {
        // do nothing as the object has no index sets
      }
    }; // class LambdaHitTestFactory
 
    template<typename Mesh_, typename Lambda_>
    MeshPart<Mesh_> make_meshpart_by_lambda_hit_test(const Mesh_& mesh, Lambda_&& lambda)
    {
      LambdaHitTestFactory<Mesh_, Lambda_> factory(mesh, std::forward<Lambda_>(lambda));
      return factory.make();
    }
 
    template<typename Mesh_, typename Lambda_>
    MeshPart<Mesh_> make_meshpart_by_filtered_lambda_hit_test(const Mesh_& mesh, const MeshPart<Mesh_>& filter, Lambda_&& lambda)
    {
      LambdaHitTestFactory<Mesh_, Lambda_> factory(mesh, filter, std::forward<Lambda_>(lambda));
      return factory.make();
    }
 
    template<typename Mesh_, typename Lambda_>
    std::unique_ptr<MeshPart<Mesh_>> make_unique_meshpart_by_lambda_hit_test(const Mesh_& mesh, Lambda_&& lambda)
    {
      LambdaHitTestFactory<Mesh_, Lambda_> factory(mesh, std::forward<Lambda_>(lambda));
      return factory.make_unique();
    }
 
    template<typename Mesh_, typename Lambda_>
    std::unique_ptr<MeshPart<Mesh_>> make_unique_meshpart_by_filtered_lambda_hit_test(const Mesh_& mesh, const MeshPart<Mesh_>& filter, Lambda_&& lambda)
    {
      LambdaHitTestFactory<Mesh_, Lambda_> factory(mesh, filter, std::forward<Lambda_>(lambda));
      return factory.make_unique();
    }
  } // namespace Geometry
} // namespace FEAT