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 <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>
 
#if defined(FEAT_HAVE_CGAL) || defined(DOXYGEN)
 
#include <cstddef>
 
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 DT_>
    class CGALWrapper
    {
    public:
      typedef DT_ DataType;
      typedef Tiny::Vector<DataType, 3> PointType;
      typedef Tiny::Matrix<DataType, 3, 3> TransformMatrix;
 
    private:
      void * _cgal_data;
 
      void _parse_mesh(std::istream & file, CGALFileMode file_mode);
 
    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;
 
      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;
 
      void transform(const TransformMatrix& trafo_mat, const PointType& translation, DataType scale = DataType(1));
 
      std::size_t bytes() const;
 
 
    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[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)