cgal.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
 
#include <iterator>
#include <kernel/base_header.hpp>
#include <kernel/geometry/mesh_part.hpp>
#include <kernel/shape.hpp>
#include <kernel/util/string.hpp>
#include <kernel/util/tiny_algebra.hpp>
#include <fstream>
 
#if defined(FEAT_HAVE_CGAL) || defined(DOXYGEN)
 
#include <cstddef>
#include <optional>
#include <functional>
 
namespace FEAT
{
  namespace Geometry
  {
 
    enum class CGALFileMode
    {
      fm_off = 0, 
      fm_obj 
    };
 
    using CGALFeature = std::vector<std::uint32_t>;
 
    using CGALFeatureNetwork = std::vector<CGALFeature>;
 
    template<typename T_>
    class CGALValueIteratorWrapper
    {
      using iterator_category = std::input_iterator_tag;
      using difference_type = Index;
      using value_type = T_;
      using pointer = T_*;
      using reference = T_&;
 
      using GeneratorType = std::function<std::optional<T_>()>;
 
    private:
      GeneratorType _generator;
      std::optional<T_> _current;
 
    public:
      explicit CGALValueIteratorWrapper(GeneratorType&& g) : _generator(std::move(g)), _current(_generator())
      {
      }
 
      T_ operator*()
      {
        return _current.value();
      }
 
      T_ operator->()
      {
        return _current.value();
      }
 
      CGALValueIteratorWrapper& operator++()
      {
        ASSERT(_current.has_value());
        _current = _generator();
        return *this;
      }
 
      CGALValueIteratorWrapper operator++(int)
      {
        ASSERT(_current.has_value());
        CGALValueIteratorWrapper tmp = *this;
        ++(*this);
        return tmp;
      }
 
      explicit operator bool() const
      {
        return static_cast<bool>(_current);
      }
 
      friend bool operator==(const CGALValueIteratorWrapper& a, const CGALValueIteratorWrapper& b)
      {
        if(!a && !b)
        {
          // Both empty
          return true;
        }
        if(!a || !b)
        {
          // One empty
          return false;
        }
        // None empty
        return true;
      }
 
      friend bool operator!=(const CGALValueIteratorWrapper& a, const CGALValueIteratorWrapper& b)
      {
        return !(a == b);
      }
    };
 
    template<typename DT_>
    class CGALWrapper;
 
    template<typename DT_>
    class CGALVerticesAroundFaceAdjactor
    {
      const CGALWrapper<DT_>* _wrapper;
    public:
 
      explicit CGALVerticesAroundFaceAdjactor(const CGALWrapper<DT_>* w) : _wrapper(w)
      {
      }
 
      Index get_num_nodes_domain() const;
 
      Index get_num_nodes_image() const;
 
      CGALValueIteratorWrapper<Index> image_begin(Index face) const;
 
      CGALValueIteratorWrapper<Index> image_end(Index face) const;
    };
 
    template<typename DT_>
    class CGALWrapper
    {
    public:
      typedef DT_ DataType;
      typedef Tiny::Vector<DataType, 3> PointType;
      typedef Tiny::Matrix<DataType, 3, 3> TransformMatrix;
 
    private:
      void * _cgal_data;
      bool _expensive_intersection = false;
 
      void _parse_mesh(std::istream & file, CGALFileMode file_mode);
 
      friend CGALVerticesAroundFaceAdjactor<DT_>;
 
    public:
      CGALWrapper(const CGALWrapper&) = delete;
      CGALWrapper& operator=(const CGALWrapper&) = delete;
      CGALWrapper(CGALWrapper&&) noexcept;
      CGALWrapper& operator=(CGALWrapper&& other) noexcept;
      virtual ~CGALWrapper();
 
      explicit CGALWrapper(const String & filename, CGALFileMode file_mode);
 
      explicit CGALWrapper(std::istream & file, CGALFileMode file_mode);
 
      explicit CGALWrapper(const std::vector<PointType>& vertices, const std::vector<std::array<Index, 3>>& faces);
 
      bool point_inside(DataType x, DataType y, DataType z) const;
 
      bool point_not_outside(DataType x, DataType y, DataType z) const;
 
      DataType squared_distance(DataType x, DataType y, DataType z) const;
 
      DataType squared_distance_to_feature(const CGALFeature& f, const PointType& p) const;
 
      PointType point(std::uint32_t idx) const;
 
      std::vector<PointType> points() const;
 
      PointType closest_point(const PointType& point) const;
 
      PointType closest_point(DataType x, DataType y, DataType z) const;
 
      PointType closest_point(const PointType& point, PointType& primitive_grad) const;
 
      PointType closest_point_on_feature(const CGALFeature& f, const PointType& p) const;
 
      CGALFeatureNetwork detect_features(DataType critical_angle);
 
      bool intersects_line(const PointType& a, const PointType& b) const;
 
      bool intersects_hexaedron(const std::array<FEAT::Geometry::CGALWrapper<DT_>::PointType, 8>& points) const;
 
      template<std::size_t num_verts_>
      bool intersects_polygon(const std::array<FEAT::Geometry::CGALWrapper<DT_>::PointType, num_verts_>& points) const
      {
        if constexpr(num_verts_ == 2)
        {
          return intersects_line(points[0], points[1]);
        }
        else if constexpr(num_verts_ == 8)
        {
          return intersects_hexaedron(points);
        }
 
        XABORTM("Not implemented");
 
        return false;
      }
 
      void transform(const TransformMatrix& trafo_mat, const PointType& translation, DataType scale = DataType(1));
 
      std::size_t bytes() const;
 
      void write_off(std::ostream& stream) const;
 
      void write_off(const String& filename) const
      {
        // try to open our output file
        const String name(filename + ".off");
        std::ofstream ofs(name.c_str());
        if(!(ofs.is_open() && ofs.good()))
        {
          throw FileError("Failed to create '" + name + "'");
        }
 
        write_off(ofs);
        ofs.close();
      }
 
      void displace_vertices(const std::vector<PointType>& offsets);
 
      Index get_num_vertices() const;
 
      Index get_num_entities(int dim) const;
 
      CGALVerticesAroundFaceAdjactor<DT_> vertices_around_face() const;
 
      std::vector<PointType> outer_normals_at_faces() const;
 
      void set_expensive_intersection(bool expensive_test);
 
    private:
      void _delete_tree();
      void _init_wrapper();
 
    }; // class CGALWrapper<typename DT_>
 
    template<typename MeshType_>
    static CGALWrapper<typename MeshType_::CoordType> cgal_wrapper_from_mesh(const MeshType_& mesh, const MeshPart<MeshType_>& part)
    {
      using MeshType = MeshType_;
      using DT = typename MeshType::CoordType;
      using ShapeType = typename MeshType::ShapeType;
      using FaceType = typename Shape::FaceTraits<ShapeType, 2>::ShapeType;
      using PointType = typename CGALWrapper<DT>::PointType;
 
      constexpr int vertices_per_face = Shape::FaceTraits<FaceType, 0>::count;
 
      static_assert(ShapeType::dimension == 3);
 
      std::vector<PointType> vertices;
      vertices.reserve(part.get_num_entities(0));
 
      // Prepare vertices
      const auto& vertex_set = mesh.get_vertex_set();
      const TargetSet& vts = part.template get_target_set<0>();
      for(Index i(0); i < part.get_num_entities(0); i++)
      {
        vertices.push_back(vertex_set[vts[i]]);
      }
 
      // Prepare faces
      std::vector<std::array<Index, 3>> faces;
 
      const auto& v_at_f = mesh.template get_index_set<2, 0>();
      const TargetSet& fts = part.template get_target_set<2>();
      const Index num_mesh_faces = part.get_num_entities(2);
 
      if constexpr(std::is_same_v<ShapeType, Shape::Hexahedron>)
      {
        // Quad faces get split into triangles
        faces.reserve(2 * num_mesh_faces);
      }
      if constexpr(std::is_same_v<ShapeType, Shape::Tetrahedron>)
      {
        // Triangle faces do not need to be split
        faces.reserve(num_mesh_faces);
      }
 
      std::array<Index, vertices_per_face> vs;
      for(Index i(0); i < num_mesh_faces; i++)
      {
        for(int j(0); j < vertices_per_face; j++)
        {
          // Get mesh index
          Index v = v_at_f(fts[i], j);
 
          // Find corresponding index in vertex target set
          for(Index k(0); k < vts.get_num_entities(); k++)
          {
            if(vts[k] == v)
            {
              vs[std::size_t(j)] = k;
              break;
            }
          }
        }
 
        if constexpr(std::is_same_v<ShapeType, Shape::Hexahedron>)
        {
          faces.push_back({vs[0], vs[1], vs[2]});
          faces.push_back({vs[3], vs[2], vs[1]});
        }
        if constexpr(std::is_same_v<ShapeType, Shape::Tetrahedron>)
        {
          faces.push_back({vs[0], vs[1], vs[2]});
        }
      }
 
      return CGALWrapper<DT>(vertices, faces);
    }
  } // namespace Geometry
} // namespace FEAT
#endif //defined(FEAT_HAVE_CGAL) || defined(DOXYGEN)