// Copyright (c) 2010, 2014 INRIA Sophia-Antipolis (France). // Copyright (c) 2017 GeometryFactory (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). // // $URL: https://github.com/CGAL/cgal/blob/v5.2/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h $ // $Id: Periodic_3_mesh_triangulation_3.h 8080601 2020-12-05T17:58:52+01:00 Mael Rouxel-Labbé // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // // Author(s) : Mikhail Bogdanov // Aymeric Pellé // Mael Rouxel-Labbé #ifndef CGAL_PERIODIC_3_MESH_3_MESH_TRIANGULATION_3_H #define CGAL_PERIODIC_3_MESH_3_MESH_TRIANGULATION_3_H #include #include // traits class #include #include #include #include // periodic triangulations #include #include // vertex and cell bases #include #include #include #include #include #include #include #include #include #include #include #include #include #include namespace CGAL { /// This class currently provides an interface between the classe /// `CGAL::Periodic_3_regular_triangulation_3` and the mesher `Mesh_3`. /// As periodic triangulations are parallelized, a lot of these functions will /// become obsolete. template class Periodic_3_regular_triangulation_3_wrapper : public Periodic_3_regular_triangulation_3 { public: typedef Sequential_tag Concurrency_tag; typedef void Lock_data_structure; void *get_lock_data_structure() const { return 0; } void set_lock_data_structure(void *) const { } typedef Periodic_3_regular_triangulation_3 Base; typedef typename Base::Geom_traits Geom_traits; typedef typename Base::Triangulation_data_structure Triangulation_data_structure; typedef typename Base::Cell_iterator Finite_cells_iterator; typedef typename Base::Facet_iterator Finite_facets_iterator; typedef typename Base::Edge_iterator Finite_edges_iterator; typedef typename Base::Vertex_iterator Finite_vertices_iterator; typedef typename Base::Vertex_handle Vertex_handle; typedef typename Base::Edge Edge; typedef typename Base::Facet Facet; typedef typename Base::Cell_handle Cell_handle; typedef typename Base::FT FT; typedef typename Base::Bare_point Bare_point; typedef typename Base::Weighted_point Weighted_point; typedef typename Base::Periodic_bare_point Periodic_bare_point; typedef typename Base::Periodic_weighted_point Periodic_weighted_point; typedef typename Base::Locate_type Locate_type; typedef typename Base::Segment Segment; typedef typename Base::Periodic_segment Periodic_segment; typedef typename Base::Triangle Triangle; typedef typename Base::Periodic_triangle Periodic_triangle; typedef typename Base::Tetrahedron Tetrahedron; typedef typename Base::Periodic_tetrahedron Periodic_tetrahedron; typedef typename Base::Offset Offset; typedef typename Base::Iso_cuboid Iso_cuboid; typedef typename Base::Conflict_tester Conflict_tester; typedef typename Base::Covering_sheets Covering_sheets; typedef typename Geom_traits::Vector_3 Vector_3; typedef typename Geom_traits::Ray_3 Ray; using Base::construct_point; using Base::construct_weighted_point; using Base::construct_segment; using Base::construct_triangle; using Base::construct_tetrahedron; using Base::domain; using Base::dual; using Base::get_offset; using Base::geom_traits; using Base::adjacent_vertices; using Base::incident_cells; using Base::incident_edges; using Base::incident_facets; using Base::is_vertex; using Base::number_of_vertices; using Base::periodic_triangle; using Base::periodic_tetrahedron; using Base::point; using Base::set_point; using Base::tds; using Base::set_offsets; #ifndef CGAL_NO_STRUCTURAL_FILTERING using Base::inexact_locate; #endif static std::string io_signature() { return Get_io_signature()(); } /// Constructor Periodic_3_regular_triangulation_3_wrapper(const Iso_cuboid& domain = Iso_cuboid(0,0,0,1,1,1), const Geom_traits& gt = Geom_traits()) : Base(domain, gt) { } /// Concurrency related template bool try_lock_cell(const Cell_handle &, int = 0) const { return true; } bool try_lock_and_get_incident_cells(Vertex_handle /* v */, std::vector& /* cells */) const { std::cerr << "ERROR: implement try_lock_and_get_incident_cells()"<< std::endl; return true; } /// Basic setters-getters // there are no infinite elements in a periodic triangulation bool is_infinite(const Vertex_handle) const { return false; } bool is_infinite(const Edge&) const { return false; } bool is_infinite(const Facet&) const { return false; } bool is_infinite(const Cell_handle) const { return false; } bool is_infinite(const Cell_handle, int) const { return false; } bool is_infinite(const Cell_handle c, int i, int j) const; Cell_handle infinite_cell() const { // there are no infinite cells in a periodic triangulation CGAL_assertion(false); return Cell_handle(); } Vertex_handle infinite_vertex() const { // there is no infinite vertex in a periodic triangulation CGAL_assertion(false); return Vertex_handle(); } int dimension() const { // there can be no degenerate dimensions in a periodic triangulation return (number_of_vertices() == 0) ? -2 : 3; } void set_domain(const Iso_cuboid& domain) { CGAL_precondition_msg(domain.xmax() - domain.xmin() == domain.ymax() - domain.ymin() && domain.xmax() - domain.xmin() == domain.zmax() - domain.zmin(), "The fundamental domain must be a cube."); Base::set_domain(domain); } Bare_point canonicalize_point(const Bare_point& p) const { return P3T3::internal::robust_canonicalize_point(p, geom_traits()); } // @todo it might be dangerous to call robust_canonicalize without also changing // = construct_periodic_point(p) (lack of consistency in the result) Weighted_point canonicalize_point(const Weighted_point& p) const { return P3T3::internal::robust_canonicalize_point(p, geom_traits()); } Triangle triangle(const Facet& f) const { Periodic_triangle ptri = periodic_triangle(f); return construct_triangle(ptri); } Tetrahedron tetrahedron(const Cell_handle c) const { Periodic_tetrahedron ptet = periodic_tetrahedron(c); return construct_tetrahedron(ptet); } /*! Copies all finite `Edge`s incident to `v` to the output iterator `edges`. Returns the resulting output iterator. Since there are no infinite edges in a periodic triangulation, this function is simply a wrapper around `incident_edges` \pre `t.dimension() > 0`, `v != Vertex_handle()`, `t.is_vertex(v)`. */ template OutputIterator finite_incident_edges(Vertex_handle v, OutputIterator edges) const { return incident_edges(v, edges); } /*! Copies the `Cell_handle`s of all finite cells incident to `v` to the output iterator `cells`. Returns the resulting output iterator. Since there are no infinite cells in a periodic triangulation, this function is simply a wrapper around `incident_cells` \pre `t.dimension() == 3`, `v != Vertex_handle()`, `t.is_vertex(v)`. */ template OutputIterator finite_incident_cells(Vertex_handle v, OutputIterator cells) const { return incident_cells(v, cells); } /*! Copies all finite `Facet`s incident to `v` to the output iterator `facets`. Returns the resulting output iterator. Since there are no infinite facets in a periodic triangulation, this function is simply a wrapper around `incident_facets` \pre `t.dimension() > 1`, `v != Vertex_handle()`, `t.is_vertex(v)`. */ template OutputIterator finite_incident_facets(Vertex_handle v, OutputIterator facets) const { return incident_facets(v, facets); } // Periodic triangulations cannot be used in parallel (yet), but the functions below // are required for compilation. Note that these functions could already // be implemented and moved to P3T3.h but to make it clear that P3M3 is // not available in parallel, they are put here and left empty. template OutputIterator incident_edges_threadsafe(Vertex_handle /* v */, OutputIterator edges) const { CGAL_assertion(false); // not yet supported return edges; } template OutputIterator incident_facets_threadsafe(Vertex_handle /* v */, OutputIterator facets) const { CGAL_assertion(false); // not yet supported return facets; } template void incident_cells_threadsafe(Vertex_handle /* v */, OutputIterator /* cells */) const { CGAL_assertion(false); // not yet supported } void clear_v_offsets() const { for (typename std::vector::iterator voit = this->v_offsets.begin(); voit != this->v_offsets.end() ; ++voit) { (*voit)->clear_offset(); } this->v_offsets.clear(); } FT compute_power_distance_to_power_sphere(const Cell_handle& c, const int i) const { typename Geom_traits::Compute_power_distance_to_power_sphere_3 cr = geom_traits().compute_power_distance_to_power_sphere_3_object(); Offset o_nb = this->neighbor_offset(c, i, c->neighbor(i)); Offset o_vt = this->get_offset(c->neighbor(i), c->neighbor(i)->index(c)); const Weighted_point& wp0 = this->point(c->vertex(0)); // need the canonical point const Weighted_point& wp1 = this->point(c->vertex(1)); const Weighted_point& wp2 = this->point(c->vertex(2)); const Weighted_point& wp3 = this->point(c->vertex(3)); const Weighted_point& wq = this->point(c->neighbor(i)->vertex(c->neighbor(i)->index(c))); const Offset& op0 = this->get_offset(c, 0); const Offset& op1 = this->get_offset(c, 1); const Offset& op2 = this->get_offset(c, 2); const Offset& op3 = this->get_offset(c, 3); const Offset& oq = o_vt - o_nb; return cr(wp0, wp1, wp2, wp3, wq, op0, op1, op2, op3, oq); } // The functions below are needed by Mesh_3 but need a specific implementation // for the periodic case because we need to try with different offsets // to get a result FT compute_power_distance_to_power_sphere(const Cell_handle& c, const Vertex_handle v) const { typename Geom_traits::Compute_power_distance_to_power_sphere_3 cr = geom_traits().compute_power_distance_to_power_sphere_3_object(); FT min_power_dist = std::numeric_limits::infinity(); for(int i = 0; i < 3; ++i) { for(int j = 0; j < 3; ++j) { for(int k = 0; k < 3; ++k) { const Weighted_point& wp0 = this->point(c->vertex(0)); // need the canonical point const Weighted_point& wp1 = this->point(c->vertex(1)); const Weighted_point& wp2 = this->point(c->vertex(2)); const Weighted_point& wp3 = this->point(c->vertex(3)); const Weighted_point& wq = this->point(v); const Offset& op0 = this->get_offset(c, 0); const Offset& op1 = this->get_offset(c, 1); const Offset& op2 = this->get_offset(c, 2); const Offset& op3 = this->get_offset(c, 3); const Offset oq(i-1, j-1, k-1); FT power_dist = cr(wp0, wp1, wp2, wp3, wq, op0, op1, op2, op3, oq); if(power_dist < min_power_dist) min_power_dist = power_dist; } } } return min_power_dist; } // Return the tetrahedron made of 'f' + 'wp' // \pre there exists an offset such that 'f.first' is in conflict with 'wp' Tetrahedron tetrahedron(const Facet& f, const Weighted_point& wp) const { Weighted_point canonic_wp = canonicalize_point(wp); Conflict_tester tester(canonic_wp, this); const Cell_handle c = f.first; const int index = f.second; for(int i = 0; i < 3; ++i) { for(int j = 0; j < 3; ++j) { for(int k = 0; k < 3; ++k) { const Offset off(i-1, j-1, k-1); if(tester(c, off)) { return construct_tetrahedron( canonic_wp, this->point(c->vertex((index+1)&3)), this->point(c->vertex((index+2)&3)), this->point(c->vertex((index+3)&3)), off, this->get_offset(c, (index+1)&3), this->get_offset(c, (index+2)&3), this->get_offset(c, (index+3)&3)); } } } } CGAL_assertion(false); return Tetrahedron(); } Bounded_side side_of_power_sphere(const Cell_handle& c, const Weighted_point& p, bool perturb = false) const { Weighted_point canonical_p = canonicalize_point(p); Bounded_side bs = ON_UNBOUNDED_SIDE; for(int i = 0; i < 3; ++i) { for(int j = 0; j < 3; ++j) { for(int k = 0; k < 3; ++k) { bs = Base::_side_of_power_sphere(c, canonical_p, Offset(i-1, j-1, k-1), perturb); if(bs != ON_UNBOUNDED_SIDE) return bs; } } } return bs; return Base::side_of_power_sphere(c, canonical_p, Offset(), perturb); } // Warning : This is a periodic version that computes the smallest possible // between 'p' and 'q', for all possible combinations of offsets FT min_squared_distance(const Bare_point& p, const Bare_point& q) const { typename Geom_traits::Compute_squared_distance_3 csd = geom_traits().compute_squared_distance_3_object(); const Bare_point cp = canonicalize_point(p); const Bare_point cq = canonicalize_point(q); FT min_sq_dist = std::numeric_limits::infinity(); for(int i = 0; i < 3; ++i) { for(int j = 0; j < 3; ++j) { for(int k = 0; k < 3; ++k) { FT sq_dist = csd(cq, construct_point(std::make_pair(cp, Offset(i-1, j-1, k-1)))); if(sq_dist < min_sq_dist) min_sq_dist = sq_dist; } } } return min_sq_dist; } // Warning : This function finds which offset 'Oq' should be applied to 'q' so // that the distance between 'p' and '(q, Oq)' is minimal. // // \pre 'p' lives in the canonical instance. Bare_point get_closest_point(const Bare_point& p, const Bare_point& q) const { Bare_point rq; const Bare_point cq = canonicalize_point(q); FT min_sq_dist = std::numeric_limits::infinity(); for(int i = -1; i < 2; ++i) { for(int j = -1; j < 2; ++j) { for(int k = -1; k < 2; ++k) { const Bare_point tcq = construct_point(std::make_pair(cq, Offset(i, j, k))); FT sq_dist = geom_traits().compute_squared_distance_3_object()(p, tcq); if(sq_dist < min_sq_dist) { rq = tcq; min_sq_dist = sq_dist; } } } } return rq; } Weighted_point get_closest_point(const Weighted_point& wp, const Weighted_point& wq) const { typename Geom_traits::Compute_weight_3 cw = geom_traits().compute_weight_3_object(); typename Geom_traits::Construct_point_3 cp = geom_traits().construct_point_3_object(); typename Geom_traits::Construct_weighted_point_3 cwp = geom_traits().construct_weighted_point_3_object(); return cwp(get_closest_point(cp(wp), cp(wq)), cw(wq)); } Triangle get_closest_triangle(const Bare_point& p, const Triangle& t) const { typename Geom_traits::Construct_vector_3 cv = geom_traits().construct_vector_3_object(); typename Geom_traits::Construct_translated_point_3 tr = geom_traits().construct_translated_point_3_object(); typename Geom_traits::Compute_squared_distance_3 csd = geom_traits().compute_squared_distance_3_object(); // It doesn't matter which point we use to canonicalize the triangle as P3M3 is necessarily // in one cover and we have to look at all the neighboring copies anyway since we do not // have control of 'p'. Bare_point canon_p0 = canonicalize_point(t[0]); Vector_3 move_to_canonical = cv(t[0], canon_p0); const std::array ct = { canon_p0, tr(t[1], move_to_canonical), tr(t[2], move_to_canonical) }; FT min_sq_dist = std::numeric_limits::infinity(); Triangle rt; for(int i = -1; i < 2; ++i) { for(int j = -1; j < 2; ++j) { for(int k = -1; k < 2; ++k) { const Triangle tt( construct_point(std::make_pair(ct[0], Offset(i, j, k))), construct_point(std::make_pair(ct[1], Offset(i, j, k))), construct_point(std::make_pair(ct[2], Offset(i, j, k)))); const FT sq_dist = csd(p, tt); if(sq_dist == FT(0)) return rt; if(sq_dist < min_sq_dist) { rt = tt; min_sq_dist = sq_dist; } } } } return rt; } // Warning: This is a periodic version that computes the smallest possible // distances between p and q, and between p and r FOR ALL POSSIBLE OFFSETS // before comparing these distances. // // It is used in facet encroachement checks in Periodic_3_mesh_3. bool greater_or_equal_power_distance(const Bare_point& p, const Weighted_point& q, const Weighted_point& r) const { typename Geom_traits::Compute_power_product_3 power_distance = geom_traits().compute_power_product_3_object(); // canonicalize the points const Weighted_point cp = geom_traits().construct_weighted_point_3_object()(canonicalize_point(p)); const Weighted_point cq = canonicalize_point(q); const Weighted_point cr = canonicalize_point(r); FT min_power_distance_to_q = std::numeric_limits::infinity(); FT min_power_distance_to_r = std::numeric_limits::infinity(); for(int i = 0; i < 3; ++i) { for(int j = 0; j < 3; ++j) { for(int k = 0; k < 3; ++k) { const Weighted_point cp_copy = construct_weighted_point(std::make_pair(cp, Offset(i-1, j-1, k-1))); const FT power_distance_to_q = power_distance(cp_copy, cq); if(power_distance_to_q < min_power_distance_to_q) min_power_distance_to_q = power_distance_to_q; const FT power_distance_to_r = power_distance(cp_copy, cr); if (power_distance_to_r < min_power_distance_to_r) min_power_distance_to_r = power_distance_to_r; } } } CGAL_postcondition(min_power_distance_to_r < 0.5 && min_power_distance_to_q < 0.5); return min_power_distance_to_q >= min_power_distance_to_r; } /// \name Locate functions /// /// Locate points within a periodic triangulation. /// /// These functions are temporarily here to interface between Mesh_3 and /// the periodic triangulations, until the latter are made parallel. /// /// \sa `CGAL::Regular_triangulation_3::locate` /// @{ bool is_vertex(const Weighted_point& p, Vertex_handle& v) const { return Base::is_vertex(canonicalize_point(p), v); } Vertex_handle nearest_power_vertex(const Bare_point& p, Cell_handle start) const { return Base::nearest_power_vertex(canonicalize_point(p), start); } /// Return the squared distance (note: _NOT_ the power distance) between the /// 'p' and the closest vertex for the power distance. std::pair nearest_power_vertex_with_sq_distance(const Bare_point& p, Cell_handle start) const { // The function below is almost a copy from 'nearest_power_vertex' in P3RT3. // Any change should be mirrored. CGAL_precondition(number_of_vertices() > 0); Bare_point canonical_p = canonicalize_point(p); Locate_type lt; int li, lj; typename Geom_traits::Compute_squared_distance_3 csd = geom_traits().compute_squared_distance_3_object(); typename Geom_traits::Construct_point_3 cp = geom_traits().construct_point_3_object(); typename Geom_traits::Construct_weighted_point_3 cwp = geom_traits().construct_weighted_point_3_object(); Offset query_offset; Cell_handle c = Base::locate(cwp(canonical_p), query_offset, lt, li, lj, start); // - start with the closest vertex from the located cell. // - repeatedly take the nearest of its incident vertices if any // - if not, we're done. // Take the opposite because periodic_locate() returns the offset such that // cell + offset contains 'p' but here we need to move 'p' query_offset = this->combine_offsets(Offset(), -query_offset); Vertex_handle nearest = Base::nearest_vertex_in_cell(c, canonical_p, query_offset); const Weighted_point& nearest_wp = this->point(nearest); Offset offset_of_nearest = Base::get_min_dist_offset(canonical_p, query_offset, nearest); FT min_sq_dist = csd(canonical_p, cp(nearest_wp), query_offset, offset_of_nearest); std::vector vs; vs.reserve(32); while(true) { Vertex_handle tmp = nearest; adjacent_vertices(nearest, std::back_inserter(vs)); for(typename std::vector::const_iterator vsit = vs.begin(); vsit != vs.end(); ++vsit) { // Can happen in 27-sheeted triangulations composed of few points if(point(*vsit) == point(nearest)) continue; const Offset min_dist_offset = this->get_min_dist_offset(canonical_p, query_offset, *vsit); if(this->compare_distance(canonical_p, point(*vsit), point(tmp), query_offset, min_dist_offset, offset_of_nearest) == SMALLER) { tmp = *vsit; offset_of_nearest = min_dist_offset; const Weighted_point& vswp = this->point(tmp); min_sq_dist = csd(canonical_p, cp(vswp), query_offset, min_dist_offset); } } if(tmp == nearest) break; vs.clear(); nearest = tmp; } return std::make_pair(nearest, min_sq_dist); } Cell_handle locate(const Weighted_point& p, Cell_handle start = Cell_handle(), bool* CGAL_assertion_code(could_lock_zone) = nullptr) const { CGAL_assertion(could_lock_zone == nullptr); return Base::locate(canonicalize_point(p), start); } Cell_handle locate(const Weighted_point& p, Vertex_handle hint, bool* CGAL_assertion_code(could_lock_zone) = nullptr) const { CGAL_assertion(could_lock_zone == nullptr); // Compared to the non-periodic version in T3, the infinite cell cannot // be used as default hint, so `Cell_handle()` is used instead. return Base::locate(canonicalize_point(p), hint == Vertex_handle() ? Cell_handle() : hint->cell()); } Cell_handle locate(const Weighted_point& p, Locate_type& l, int& i, int& j, Cell_handle start = Cell_handle(), bool* CGAL_assertion_code(could_lock_zone) = nullptr) const { CGAL_assertion(could_lock_zone == nullptr); return Base::locate(canonicalize_point(p), l, i, j, start); } Cell_handle locate(const Weighted_point& p, Locate_type& l, int& i, int& j, Vertex_handle hint, bool* CGAL_assertion_code(could_lock_zone) = nullptr) const { CGAL_assertion(could_lock_zone == nullptr); return Base::locate(canonicalize_point(p), l, i, j, hint == Vertex_handle() ? Cell_handle() : hint->cell()); } /// @} /// \name Conflict functions /// Returns the vertices on the interior of the conflict hole. /// /// These functions are temporarily here to interface between Mesh_3 and /// the periodic triangulations, until the latter are made parallel. /// /// @{ template OutputIterator vertices_inside_conflict_zone(const Weighted_point& /* p */, Cell_handle /* c */, OutputIterator res) const { return res; CGAL_assertion(false); // not yet supported } template Triple find_conflicts(const Weighted_point& p, Cell_handle c, OutputIteratorBoundaryFacets bfit, OutputIteratorCells cit, OutputIteratorInternalFacets ifit, bool* CGAL_assertion_code(could_lock_zone) = nullptr, const Facet* /* this_facet_must_be_in_the_cz */ = nullptr, bool* /* the_facet_is_in_its_cz */ = nullptr) const { CGAL_triangulation_precondition(could_lock_zone == nullptr); CGAL_triangulation_precondition(number_of_vertices() != 0); clear_v_offsets(); Weighted_point canonic_p = canonicalize_point(p); Locate_type lt; int li, lj; c = locate(canonic_p, lt, li, lj, c); std::vector facets; facets.reserve(64); std::vector cells; cells.reserve(32); Conflict_tester tester(canonic_p, this); Triple >, typename std::back_insert_iterator >, OutputIteratorInternalFacets> tit = Base::find_conflicts(c, tester, make_triple(std::back_inserter(facets), std::back_inserter(cells), ifit)); ifit = tit.third; // Reset the conflict flag on the boundary. for(typename std::vector::iterator fit=facets.begin(); fit != facets.end(); ++fit) { fit->first->neighbor(fit->second)->tds_data().clear(); *bfit++ = *fit; } // Reset the conflict flag in the conflict cells. for(typename std::vector::iterator ccit=cells.begin(); ccit != cells.end(); ++ccit) { (*ccit)->tds_data().clear(); *cit++ = *ccit; } return make_triple(bfit, cit, ifit); } template std::pair find_conflicts(const Weighted_point &p, Cell_handle c, OutputIteratorBoundaryFacets bfit, OutputIteratorCells cit, bool* could_lock_zone = nullptr) const { CGAL_assertion(could_lock_zone == nullptr); Triple t = find_conflicts(p, c, bfit, cit, Emptyset_iterator(), could_lock_zone); return std::make_pair(t.first, t.second); } /// @} void set_point(const Vertex_handle v, const Vector_3& move, const Weighted_point& new_position) { // calling robust canonical here means we don't necessarily have // canonical(v + move) = new_position... @fixme return Base::set_point(v, move, canonicalize_point(new_position)); } /// \name Insert functions /// /// Insert points in the triangulation. /// /// These functions are temporarily here to interface between Mesh_3 and /// the periodic triangulations, until the latter are made parallel. /// /// @{ template Vertex_handle insert_in_hole(const Weighted_point& p, CellIt /*cell_begin*/, CellIt /*cell_end*/, Cell_handle /*begin*/, int /*i*/) { // Could be optimized to use the conflict zone by extracting relevant // code from periodic_insert() in P3T3. return this->insert(p); } Vertex_handle insert(const Weighted_point& p, Cell_handle start = Cell_handle(), bool* CGAL_assertion_code(could_lock_zone) = nullptr) { CGAL_assertion(could_lock_zone == nullptr); return Base::insert(canonicalize_point(p), start); } Vertex_handle insert(const Weighted_point& p, Vertex_handle hint, bool* CGAL_assertion_code(could_lock_zone) = nullptr) { CGAL_assertion(could_lock_zone == nullptr); // compared to the non-periodic version in T3, the infinite cell cannot // be used; `Cell_handle()` is used instead return Base::insert(canonicalize_point(p), hint == Vertex_handle() ? Cell_handle() : hint->cell()); } Vertex_handle insert(const Weighted_point& p, Locate_type lt, Cell_handle loc, int li, int lj, bool* CGAL_assertion_code(could_lock_zone) = nullptr) { CGAL_assertion(could_lock_zone == nullptr); return Base::insert(canonicalize_point(p), lt, loc, li, lj); } /// @} /// Remove functions bool remove(Vertex_handle v, bool* CGAL_assertion_code(could_lock_zone) = nullptr) { CGAL_assertion(could_lock_zone == nullptr); bool b = Base::remove_if_no_cover_change(v); CGAL_postcondition(this->is_1_cover()); // do not ever allow cover change return b; } /// Dual computations Object dual(const Facet & f) const { Segment s = construct_segment(Base::dual(f)); return make_object(s); } void dual_segment(const Facet& facet, Bare_point& p, Bare_point& q) const { typename Base::Periodic_segment_3 ps = Base::dual(facet); p = construct_point(ps[0]); q = construct_point(ps[1]); } void dual_segment_exact(const Facet& facet, Bare_point& p, Bare_point& q) const { typedef typename Kernel_traits::Kernel Kernel; typedef Exact_predicates_exact_constructions_kernel EKernel; typedef Cartesian_converter To_exact; typedef Cartesian_converter Back_from_exact; typedef CGAL::Periodic_3_regular_triangulation_traits_3 Exact_Rt; typedef typename Exact_Rt::Point_3 EPoint_3; To_exact to_exact; Back_from_exact back_from_exact; Exact_Rt etraits(to_exact(domain())); const typename EKernel::Iso_cuboid_3& dom = etraits.get_domain(); Exact_Rt::Construct_weighted_circumcenter_3 exact_weighted_circumcenter = etraits.construct_weighted_circumcenter_3_object(); Exact_Rt::Construct_point_3 exact_construct_point = etraits.construct_point_3_object(); Cell_handle c = facet.first; const int i = facet.second; Cell_handle n = c->neighbor(i); EPoint_3 exact_wc1 = exact_weighted_circumcenter( to_exact(point(c->vertex(0))), to_exact(point(c->vertex(1))), to_exact(point(c->vertex(2))), to_exact(point(c->vertex(3))), get_offset(c, 0), get_offset(c, 1), get_offset(c, 2), get_offset(c, 3)); EPoint_3 exact_wc2 = exact_weighted_circumcenter( to_exact(point(n->vertex(0))), to_exact(point(n->vertex(1))), to_exact(point(n->vertex(2))), to_exact(point(n->vertex(3))), get_offset(n, 0), get_offset(n, 1), get_offset(n, 2), get_offset(n, 3)); typename EKernel::Point_3 dp; // get the offset of the first weighted circumcenter Offset transl_wc1; while(true) /* while not in */ { dp = etraits.construct_point_3_object()(exact_wc1, transl_wc1); if(dp.x() < dom.xmin()) transl_wc1.x() += 1; else if(dp.y() < dom.ymin()) transl_wc1.y() += 1; else if(dp.z() < dom.zmin()) transl_wc1.z() += 1; else if(!(dp.x() < dom.xmax())) transl_wc1.x() -= 1; else if(!(dp.y() < dom.ymax())) transl_wc1.y() -= 1; else if(!(dp.z() < dom.zmax())) transl_wc1.z() -= 1; else break; } // get the offset of the second weighted circumcenter Offset transl_wc2; while(true) /* while not in */ { dp = etraits.construct_point_3_object()(exact_wc2, transl_wc2); if(dp.x() < dom.xmin()) transl_wc2.x() += 1; else if(dp.y() < dom.ymin()) transl_wc2.y() += 1; else if(dp.z() < dom.zmin()) transl_wc2.z() += 1; else if(!(dp.x() < dom.xmax())) transl_wc2.x() -= 1; else if(!(dp.y() < dom.ymax())) transl_wc2.y() -= 1; else if(!(dp.z() < dom.zmax())) transl_wc2.z() -= 1; else break; } Offset off = this->neighbor_offset(c, i, n); Offset o1 = -transl_wc1; Offset o2 = this->combine_offsets(-transl_wc2, -off); Offset cumm_off((std::min)(o1.x(), o2.x()), (std::min)(o1.y(), o2.y()), (std::min)(o1.z(), o2.z())); EPoint_3 ewc1 = exact_construct_point(exact_wc1, transl_wc1); EPoint_3 ewc2 = exact_construct_point(exact_wc2, transl_wc2); p = back_from_exact(exact_construct_point(ewc1, o1 - cumm_off)); q = back_from_exact(exact_construct_point(ewc2, o2 - cumm_off)); } // dual rays are impossible in a periodic triangulation since there are no // infinite cells, but these functions are still required to compile Mesh_3 void dual_ray(const Facet& /*f*/, Ray& /*ray*/) const { CGAL_assertion(false); } void dual_ray_exact(const Facet& /*facet*/, Ray& /*ray*/) const { CGAL_assertion(false); } }; namespace details { template struct Periodic_3_mesh_geom_traits_generator { private: // - 'Robust_periodic_weighted_circumcenter_traits_3': provides robust periodic // constructions for construct_weighted_circumcenter() // - 'Periodic_3_regular_triangulation_traits_3': makes the base (non-periodic) // traits usable in a periodic setting // - 'Robust_weighted_circumcenter_filtered_traits_3': provides robust versions // of some non-periodic constructions. typedef Robust_periodic_weighted_circumcenter_traits_3< Periodic_3_regular_triangulation_traits_3< Robust_weighted_circumcenter_filtered_traits_3 > > Geom_traits; public: typedef Geom_traits type; typedef type Type; }; // end struct Periodic_3_mesh_geom_traits_generator } // end namespace details template class Periodic_3_mesh_triangulation_3 { // default K typedef typename Default::Get::Kernel>::type K; // traits typedef typename details::Periodic_3_mesh_geom_traits_generator::type Geom_traits; // Periodic vertex and cell bases typedef Periodic_3_triangulation_ds_vertex_base_3<> VbDS; typedef Regular_triangulation_vertex_base_3 PVb; typedef Periodic_3_triangulation_ds_cell_base_3<> CbDS; typedef Regular_triangulation_cell_base_3 RCb; typedef Regular_triangulation_cell_base_with_weighted_circumcenter_3 PCb; typedef Mesh_vertex_base_3 Default_Vb; typedef Mesh_cell_base_3 Default_Cb; // default Vb/Cb typedef typename Default::Get::type Vertex_base; typedef typename Default::Get::type Cell_base; // Triangulation and tds typedef Triangulation_data_structure_3 Tds; typedef Periodic_3_regular_triangulation_3_wrapper Triangulation; public: typedef Triangulation type; typedef type Type; }; } // namespace CGAL #endif // CGAL_PERIODIC_3_MESH_3_MESH_TRIANGULATION_3_H