diff --git a/CGAL3DPolyhedralMesher/CMakeLists.txt b/CGAL3DPolyhedralMesher/CMakeLists.txt new file mode 100644 index 0000000000..53ddae4de6 --- /dev/null +++ b/CGAL3DPolyhedralMesher/CMakeLists.txt @@ -0,0 +1,46 @@ +if (NOT ParaView_BINARY_DIR) + find_package(ParaView REQUIRED) +endif() + +#Find CGAL +find_package(CGAL) +if(CGAL_FOUND) + include( ${CGAL_USE_FILE} ) +endif(CGAL_FOUND) + +if(ParaView_FOUND) + #If the plugin is used internally we don't need to include. + if (NOT ParaView_BINARY_DIR) + include(${PARAVIEW_USE_FILE}) + endif(NOT ParaView_BINARY_DIR) +else(ParaView_FOUND) + if(VTK_INCLUDE_DIRS) + INCLUDE_DIRECTORIES(${VTK_INCLUDE_DIRS}) + endif(VTK_INCLUDE_DIRS) +endif(ParaView_FOUND) + +#Find CGAL +find_package(CGAL) +if(CGAL_FOUND) + include( ${CGAL_USE_FILE} ) +endif(CGAL_FOUND) + +if(ParaView_FOUND OR VTK_INCLUDE_DIRS) + message(NOTICE "inside") +# create a paraview plugin containing server manager xml and the server +# manager classes to build +# this plugin can be loaded on the server side + + #Paraview's macro to add the plugin. It takes care of all the vtk + #and paraview parts of the process, like link and integration + #in the UI + ADD_PARAVIEW_PLUGIN(CGAL3DPolyhedralMesher "1.0" + SERVER_MANAGER_XML vtkCGAL3DPolyhedralMesher.xml + SERVER_MANAGER_SOURCES vtkCGAL3DPolyhedralMesher.cxx) +endif(ParaView_FOUND OR VTK_INCLUDE_DIRS) + +if(CGAL_FOUND) + target_link_libraries(CGAL3DPolyhedralMesher PRIVATE + CGAL::CGAL + ${Boost_LIBRARIES}) +endif (CGAL_FOUND) \ No newline at end of file diff --git a/CGAL3DPolyhedralMesher/plugin.cmake b/CGAL3DPolyhedralMesher/plugin.cmake new file mode 100644 index 0000000000..80f07c745f --- /dev/null +++ b/CGAL3DPolyhedralMesher/plugin.cmake @@ -0,0 +1,4 @@ +pv_plugin(CGAL3DPolyhedralMesher + DESCRIPTION "CGAL 3D Polyhedral Mesher" + AUTOLOAD + DEFAULT_ENABLED) diff --git a/CGAL3DPolyhedralMesher/vtkCGAL3DPolyhedralMesher.cxx b/CGAL3DPolyhedralMesher/vtkCGAL3DPolyhedralMesher.cxx new file mode 100644 index 0000000000..91a24b98e0 --- /dev/null +++ b/CGAL3DPolyhedralMesher/vtkCGAL3DPolyhedralMesher.cxx @@ -0,0 +1,313 @@ +/** +* \class vtkCGAL3DPolyhedralMesher +* +* \brief Generates a polyhedral mesh of the bounding domain based on features related to the interior surfaces. +* ODT and Lloyd optimization methods are supported. Exude and perturbe options are also supported, but can +* have an impact on the refinement criteria. See the official CGAL documentation for more information. +* +* Inputs: Interior Surfaces (port 0, vtkPolyData), Bounding Domain (port 1, vtkPolyData) +* Output: Polyhedral Domain with features (port 0, vtkUnstructuredGrid) +*/ + +#include "vtkCGAL3DPolyhedralMesher.h" + +// VTK +#include +#include +#include +#include +#include +#include + +#include + +// CGAL +#include +#include +#include +#include +#include +#include +#include +#include +#include + +vtkStandardNewMacro(vtkCGAL3DPolyhedralMesher); + +// ---------------------------------------------------------------------------- +vtkCGAL3DPolyhedralMesher::vtkCGAL3DPolyhedralMesher() +{ + this->SetNumberOfInputPorts(2); + this->SetNumberOfOutputPorts(1); + + // Mesh Criteria + this->EdgeSize = 0.025; + this->FacetAngle = 25; + this->FacetSize = 0.05; + this->FacetDistance = 0.005; + this->CellRadiusEdgeRatio = 3; + this->CellSize = 0.05; + + // Optimizer flags + this->Lloyd = false; + this->Odt = false; + this->Perturb = false; + this->Exude = false; + + this->TopologicalStructure = vtkCGAL3DPolyhedralMesher::TopologicalStructures::MANIFOLD; + this->ConfineSurfacePoints = true; +} + +//---------------------------------------------------------------------------- +void vtkCGAL3DPolyhedralMesher::PrintSelf(ostream& os, vtkIndent indent) +{ + // TODO: Write PrintSelf + this->Superclass::PrintSelf(os, indent); +} + +//---------------------------------------------------------------------------- + +/** @brief Getter for interior surface +* +* @return vtkPolyData* The interior surface features +*/ +vtkPolyData* vtkCGAL3DPolyhedralMesher::GetInteriorSurfaces() +{ + if (this->GetNumberOfInputConnections(0) < 1) { + return nullptr; + } + + return vtkPolyData::SafeDownCast(this->GetInputDataObject(0, 0)); +} + +//---------------------------------------------------------------------------- + +/** @brief Getter for Bounding Domain +* +* @return vtkPolyData* +*/ +vtkPolyData* vtkCGAL3DPolyhedralMesher::GetBoundingDomain() +{ + if (this->GetNumberOfInputConnections(1) < 1) { + return nullptr; + } + + return vtkPolyData::SafeDownCast(this->GetInputDataObject(1, 0)); +} + +//---------------------------------------------------------------------------- + +int vtkCGAL3DPolyhedralMesher::RequestData(vtkInformation* vtkNotUsed(request), + vtkInformationVector** inputVector, + vtkInformationVector* outputVector) +{ + vtkPolyData* inputInteriorSurfaces = vtkPolyData::GetData(inputVector[0]->GetInformationObject(0)); + vtkPolyData* inputBoundingDomain = vtkPolyData::GetData(inputVector[1]->GetInformationObject(0)); + vtkUnstructuredGrid* output = vtkUnstructuredGrid::GetData(outputVector->GetInformationObject(0)); + + if (inputInteriorSurfaces == nullptr) + { + vtkErrorMacro("Interior Surfaces input is null."); + return 0; + } + + if (inputInteriorSurfaces->GetNumberOfPoints() == 0) + { + vtkErrorMacro("Interior Surfaces input does not contain any points."); + return 0; + } + + if (inputInteriorSurfaces->GetPolys() == nullptr) + { + vtkErrorMacro("Interior Surfaces input does not contain any cell structure."); + return 0; + } + + if (inputBoundingDomain == nullptr) + { + vtkErrorMacro("Bounding Domain input is null."); + return 0; + } + + if (inputBoundingDomain->GetNumberOfPoints() == 0) + { + vtkErrorMacro("Bounding Domain input does not contain any points."); + return 0; + } + + if (inputBoundingDomain->GetPolys() == nullptr) + { + vtkErrorMacro("Bounding Domain input does not cell structure."); + return 0; + } + + // Start + this->SetProgressText("CGAL Polyhedral Mesher"); + this->UpdateProgress(0); + + // Convert from VTK to CGAL + vtkTimerLog::MarkStartEvent("VTK To CGAL conversion"); + + typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + typedef CGAL::Mesh_polyhedron_3::type Polyhedron; + Polyhedron* interiorSurfaces = new Polyhedron(); + Polyhedron* boundingDomain = new Polyhedron(); + this->vtkPolyDataToPolygonMesh(inputInteriorSurfaces, *interiorSurfaces); + this->vtkPolyDataToPolygonMesh(inputBoundingDomain, *boundingDomain); + + vtkTimerLog::MarkEndEvent("VTK To CGAL conversion"); + this->UpdateProgress(0.1); + + // Generate mesh + vtkTimerLog::MarkStartEvent("Mesh Generation"); + + // Domain + typedef CGAL::Polyhedral_mesh_domain_with_features_3 Mesh_domain; + Mesh_domain domain(*interiorSurfaces, *boundingDomain); + + // Criteria + typedef CGAL::Mesh_triangulation_3::type Tr; + typedef CGAL::Mesh_criteria_3 Mesh_criteria; + Mesh_criteria criteria( CGAL::parameters::edge_size = this->EdgeSize, + CGAL::parameters::facet_angle = this->FacetAngle, + CGAL::parameters::facet_size = this->FacetSize, + CGAL::parameters::facet_distance = this->FacetDistance, + CGAL::parameters::cell_radius_edge_ratio = this->CellRadiusEdgeRatio, + CGAL::parameters::cell_size = this->CellSize); + + // C3t3 init + typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3; + C3t3 c3t3; + + if (this->ConfineSurfacePoints) + { + CGAL::Mesh_3::internal::init_c3t3_with_features(c3t3, domain, criteria); + typedef C3t3::Triangulation::Point Weighted_point; + double pt[3] = { 0.0, 0.0, 0.0 }; + for (vtkIdType i = 0; i < inputInteriorSurfaces->GetNumberOfPoints(); ++i) + { + inputInteriorSurfaces->GetPoint(i, pt); + const Weighted_point p(Weighted_point::Point(pt[0], pt[1], pt[2])); + Tr::Vertex_handle vh = c3t3.triangulation().insert(p); + c3t3.add_to_complex(vh, 0); + } + } + else + { + domain.detect_features(); + c3t3 = CGAL::make_mesh_3(domain, criteria); + } + + CGAL::parameters::internal::Manifold_options manifoldOption; + switch (this->TopologicalStructure) { + case vtkCGAL3DPolyhedralMesher::TopologicalStructures::MANIFOLD: + { + manifoldOption = CGAL::parameters::manifold(); + break; + } + case vtkCGAL3DPolyhedralMesher::TopologicalStructures::MANIFOLD_WITH_BOUNDARY: + { + manifoldOption = CGAL::parameters::manifold_with_boundary(); + break; + } + case vtkCGAL3DPolyhedralMesher::TopologicalStructures::NON_MANIFOLD: + { + manifoldOption = CGAL::parameters::non_manifold(); + break; + } + default: + { + vtkErrorMacro("Unknown option."); + return 0; + } + } + + CGAL::refine_mesh_3(c3t3, domain, criteria, + CGAL::parameters::internal::Lloyd_options(this->Lloyd), + CGAL::parameters::internal::Odt_options(this->Odt), + CGAL::parameters::internal::Perturb_options(this->Perturb), + CGAL::parameters::internal::Exude_options(this->Exude), + CGAL::parameters::internal::Manifold_options(manifoldOption) // or manifold_with_boundary() or non_manifold() + ); + + vtkTimerLog::MarkEndEvent("Mesh Generation"); + + // Convert back to VTK + vtkTimerLog::MarkStartEvent("CGAL to VTK conversion"); + vtkNew result; + CGAL::output_c3t3_to_vtk_unstructured_grid(c3t3, result); + output->ShallowCopy(result); + vtkTimerLog::MarkEndEvent("CGAL to VTK conversion"); + + // Cleaning + delete interiorSurfaces; + delete boundingDomain; + interiorSurfaces = nullptr; + boundingDomain = nullptr; + + return 1; +} + +// ------------------------------------------------------------------------------ +int vtkCGAL3DPolyhedralMesher::FillOutputPortInformation(int, vtkInformation *info) +{ + info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkUnstructuredGrid"); + return 1; +} + +//---------------------------------------------------------------------------- + +/** @brief Converts a vtkPolyData (VTK) into a Polygonal Mesh (CGAL). Code taken from the +* Polyhedron demo located at https://github.com/CGAL/cgal/blob/master/Polyhedron/demo/Polyhedron/Plugins/IO/VTK_io_plugin.cpp +* This method does not write into our PolyData structure. Hence, we do not need to copy them. +* +* @param polyData The input PolyData +* @param tmesh The resulting Polygon Mesh +* @return bool Success (true) or failure (false) +*/ +template +bool vtkCGAL3DPolyhedralMesher::vtkPolyDataToPolygonMesh(vtkPointSet* polyData, TM& tmesh) +{ + typedef typename boost::property_map::type VPMap; + typedef typename boost::property_map_value::type Point_3; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + VPMap vpmap = get(CGAL::vertex_point, tmesh); + + // get nb of points and cells + vtkIdType nb_points = polyData->GetNumberOfPoints(); + vtkIdType nb_cells = polyData->GetNumberOfCells(); + + //extract points + std::vector vertex_map(nb_points); + for (vtkIdType i = 0; i < nb_points; ++i) + { + double coords[3]; + polyData->GetPoint(i, coords); + + vertex_descriptor v = add_vertex(tmesh); + put(vpmap, v, Point_3(coords[0], coords[1], coords[2])); + vertex_map[i] = v; + } + + //extract cells + for (vtkIdType i = 0; i < nb_cells; ++i) + { + if (polyData->GetCellType(i) != 5 + && polyData->GetCellType(i) != 7 + && polyData->GetCellType(i) != 9) //only supported cells are triangles, quads and polygons + continue; + vtkCell* cell_ptr = polyData->GetCell(i); + + vtkIdType nb_vertices = cell_ptr->GetNumberOfPoints(); + if (nb_vertices < 3) + return false; + std::vector vr(nb_vertices); + for (vtkIdType k = 0; k < nb_vertices; ++k) + vr[k] = vertex_map[cell_ptr->GetPointId(k)]; + + CGAL::Euler::add_face(vr, tmesh); + } + + return true; +} \ No newline at end of file diff --git a/CGAL3DPolyhedralMesher/vtkCGAL3DPolyhedralMesher.h b/CGAL3DPolyhedralMesher/vtkCGAL3DPolyhedralMesher.h new file mode 100644 index 0000000000..a00c585937 --- /dev/null +++ b/CGAL3DPolyhedralMesher/vtkCGAL3DPolyhedralMesher.h @@ -0,0 +1,179 @@ +#ifndef vtkCGAL3DPolyhedralMesher_h +#define vtkCGAL3DPolyhedralMesher_h + +#include "vtkPolyDataAlgorithm.h" + +class vtkPointSet; + +class vtkCGAL3DPolyhedralMesher : public vtkPolyDataAlgorithm +{ +public: + static vtkCGAL3DPolyhedralMesher* New(); + vtkTypeMacro(vtkCGAL3DPolyhedralMesher, vtkPolyDataAlgorithm); + void PrintSelf(ostream& os, vtkIndent indent) override; + + //@{ + /** + * This constant is used as an upper bound for the distance + * between two protecting ball centers that are consecutive on a 1-feature. + * This parameter has to be set to a positive value when 1-dimensional features protection is used. + */ + vtkSetMacro(EdgeSize, double); + vtkGetMacro(EdgeSize, double); + //@} + + //@{ + /** + * This parameter controls the shape of surface facets. Specifically, it is a lower bound for the angle (in degrees) + * of surface facets. When boundary surfaces are smooth, the termination of the meshing process is guaranteed + * if this angular bound is at most 30 degrees. + */ + vtkSetMacro(FacetAngle, double); + vtkGetMacro(FacetAngle, double); + //@} + + //@{ + /** + * This parameter controls the size of surface facets. + * Each surface facet has a surface Delaunay ball which is a ball circumscribing the surface facet + * and centered on the surface patch. The parameter facet_size is providing an upper bound + * for the radii of surface Delaunay balls. + */ + vtkSetMacro(FacetSize, double); + vtkGetMacro(FacetSize, double); + //@} + + //@{ + /** + * This parameter controls the approximation error of boundary and subdivision surfaces. + * It provides an upper bound for the distance between the circumcenter of a surface facet + * and the center of a surface Delaunay ball of this facet. + */ + vtkSetMacro(FacetDistance, double); + vtkGetMacro(FacetDistance, double); + //@} + + //@{ + /** + * This parameter controls the shape of mesh cells (but can't filter slivers, as we discussed earlier). + * It is an upper bound for the ratio between the circumradius of a mesh tetrahedron and its shortest edge. + * There is a theoretical bound for this parameter: the Delaunay refinement process is guaranteed to terminate + * for values of cell_radius_edge_ratio bigger than 2. + */ + vtkSetMacro(CellRadiusEdgeRatio, double); + vtkGetMacro(CellRadiusEdgeRatio, double); + //@} + + //@{ + /** + * This scalar parameter controls the size of mesh tetrahedra. + * It provides an upper bound on the circumradii of the mesh tetrahedra. + */ + vtkSetMacro(CellSize, double); + vtkGetMacro(CellSize, double); + //@} + + //@{ + /** + * Activates/deactivates the Lloyd smoother. + */ + vtkSetMacro(Lloyd, bool); + vtkGetMacro(Lloyd, bool); + vtkBooleanMacro(Lloyd, bool); + //@} + + //@{ + /** + * Activates/deactivates the ODT smoother. + */ + vtkSetMacro(Odt, bool); + vtkGetMacro(Odt, bool); + vtkBooleanMacro(Odt, bool); + //@} + + //@{ + /** + * Activates/deactivates the perturber. + */ + vtkSetMacro(Perturb, bool); + vtkGetMacro(Perturb, bool); + vtkBooleanMacro(Perturb, bool); + //@} + + //@{ + /** + * Activates/deactivates the exuder. + */ + vtkSetMacro(Exude, bool); + vtkGetMacro(Exude, bool); + vtkBooleanMacro(Exude, bool); + //@} + + typedef enum { + MANIFOLD = 1, + MANIFOLD_WITH_BOUNDARY, + NON_MANIFOLD + } TopologicalStructures; + + //@{ + /** + * Specifies the manifold of the interior surfaces + */ + vtkGetMacro(TopologicalStructure, int); + vtkSetMacro(TopologicalStructure, int); + //@} + + virtual void SetTopologicalStructureToManifold(void) + { this->SetTopologicalStructure(vtkCGAL3DPolyhedralMesher::TopologicalStructures::MANIFOLD); } + + virtual void SetTopologicalStructureToManifoldWithBoundary(void) + { this->SetTopologicalStructure(vtkCGAL3DPolyhedralMesher::TopologicalStructures::MANIFOLD_WITH_BOUNDARY); } + + virtual void SetTopologicalStructureToNonManifold(void) + { this->SetTopologicalStructure(vtkCGAL3DPolyhedralMesher::TopologicalStructures::NON_MANIFOLD); } + + //@{ + /** + * If true, every point of the interior surface will be part of the tetrahedral mesh output. + * If false, the default behavior is applied. + */ + vtkSetMacro(ConfineSurfacePoints, bool); + vtkGetMacro(ConfineSurfacePoints, bool); + vtkBooleanMacro(ConfineSurfacePoints, bool); + //@} + + vtkPolyData* GetInteriorSurfaces(); + vtkPolyData* GetBoundingDomain(); + +protected: + vtkCGAL3DPolyhedralMesher(); + ~vtkCGAL3DPolyhedralMesher(){} + + virtual int RequestData(vtkInformation*, vtkInformationVector**, vtkInformationVector*) override; + virtual int FillOutputPortInformation(int, vtkInformation *info) override; + +private: + vtkCGAL3DPolyhedralMesher(const vtkCGAL3DPolyhedralMesher&) = delete; + void operator=(const vtkCGAL3DPolyhedralMesher&) = delete; + + template + bool vtkPolyDataToPolygonMesh(vtkPointSet* poly_data, TM& tmesh); + + // Mesh Criteria + double EdgeSize; + double FacetAngle; + double FacetSize; + double FacetDistance; + double CellRadiusEdgeRatio; + double CellSize; + + // Optimizer flags + bool Lloyd; + bool Odt; + bool Perturb; + bool Exude; + + int TopologicalStructure; + bool ConfineSurfacePoints; +}; +#endif diff --git a/CGAL3DPolyhedralMesher/vtkCGAL3DPolyhedralMesher.xml b/CGAL3DPolyhedralMesher/vtkCGAL3DPolyhedralMesher.xml new file mode 100644 index 0000000000..ae9c56e7e4 --- /dev/null +++ b/CGAL3DPolyhedralMesher/vtkCGAL3DPolyhedralMesher.xml @@ -0,0 +1,213 @@ + + + + + + + + + + + Interior surface part of the bounding domain. + + + + + + + + + + + + + Bounding box defining the polyhedral domain. + + + + + + + + + + + + + This constant is used as an upper bound for the distance between two protecting + ball centers that are consecutive on a 1-feature. This parameter has to be set + to a positive value when 1-dimensional features protection is used. + + + + + + This parameter controls the shape of surface facets. Specifically, it is a lower bound + for the angle (in degrees) of surface facets. When boundary surfaces are smooth, + the termination of the meshing process is guaranteed if this angular bound is at most 30 degrees. + + + + + + This parameter controls the size of surface facets. Each surface facet has a surface Delaunay ball + which is a ball circumscribing the surface facet and centered on the surface patch. + The parameter facet_size is providing an upper bound for the radii of surface Delaunay balls. + + + + + + This parameter controls the approximation error of boundary and subdivision surfaces. It provides + an upper bound for the distance between the circumcenter of a surface facet and the center of a + surface Delaunay ball of this facet. + + + + + + This parameter controls the shape of mesh cells (but can't filter slivers, as we discussed earlier). + It is an upper bound for the ratio between the circumradius of a mesh tetrahedron and its shortest edge. + There is a theoretical bound for this parameter: the Delaunay refinement process is guaranteed to terminate + for values of cell_radius_edge_ratio bigger than 2. + + + + + + This scalar parameter controls the size of mesh tetrahedra. + It provides an upper bound on the circumradii of the mesh tetrahedra. + + + + + + + + + + + + + + + + Activates/deactivates the Lloyd smoother. + + + + + + + Activates/deactivates the ODT smoother. + + + + + + + Activates/deactivates the perturber. + + + + + + + Activates/deactivates the exuder. + + + + + + + + + + + + + + + + + + Specifies the manifold of the interior surfaces + + + + + + + If true, every point of the interior surface will be part of the tetrahedral mesh output. + If false, the default behavior is applied. + + + + + + + + + +