#ifndef OERSTED_MESH_H #define OERSTED_MESH_H #include #include #include #include #include "Eigen" #include "Eigen/Sparse" #include "Eigen/SparseLU" #include "Sketch.hpp" #include "BoundaryConstraint.h" #include "MappedBoundaryPair.h" #include "Edge.h" #include "InsertionQueuer.h" #include "Point.h" #include "DartTriangle.h" enum class LocateTriangleResult { Interior, Exterior, Boundary, Point // TODO: Enumerate cases in function when triangle is located by a point near the boundary }; enum class InsertPointResult { Success, Midpoint, Duplicate, Failure }; class Mesh { // TODO: Namespaces. Also, need to rename a bunch of things to be more consistent and clear w.r.t. how various arrays are accessed (e.g. by edge/dart index versus plain index) public: Mesh() {}; Mesh(Sketch &s); friend class BoundaryConstraint; friend class MappedBoundaryPair; friend class CircumcenterQueuer; friend class MidpointQueuer; double MinimumElementQuality = 0.0; double MinimumElementSize = 0.0; double MaximumElementSize = DBL_MAX; bool are_intersecting(size_t ei, size_t ej) const; bool edges_are_optimal() const; bool edges_are_valid() const; bool in_triangle(Point const p, size_t ei) const; bool is_constrained(size_t ei) const { return Edges[ei].Constraint != 0; }; bool is_encroached(Point const p, size_t ei) const; bool is_locally_optimal(size_t ei) const; bool is_protruding(size_t ei) const; bool is_valid(size_t ei) const; bool orientation(size_t ei) const { return Edges[ei].Orientation; }; bool refine(); bool refine_once(); double circumradius(size_t ei) const; double length(size_t ei) const; double shortest_edge_length(size_t ei) const; size_t next(size_t ei) const { return Edges[ei].Next; }; size_t node(size_t ei) const { return Edges[ei].Node; }; size_t node(Edge const e) const { return e.Node; }; size_t num_points() const { return Points.size(); }; size_t num_edges() const; size_t num_triangles() const { return Triangles.size(); }; size_t prev(size_t ei) const { return Edges[ei].Prev; }; size_t size_dart_constraints() const {return DartConstraints.size(); }; size_t size_edges() const { return Edges.size(); }; size_t size_points() const { return Points.size(); }; size_t size_triangles() const { return Triangles.size(); }; size_t size_contours() const { return Contours.size(); }; size_t size_boundary_constraints() const { return BoundaryConstraints.size(); }; size_t twin(size_t ei) const { return Edges[ei].Twin; }; void add_to_queue(std::unique_ptr ic) { Queue.push_back(std::move(ic)); }; void create(); void save_as(std::string path, std::string file_name) const; void smooth(); // TODO: rename these methods something like, dart_constraint_from_edge_index, curve_from_edge_index DartConstraint const dart_constraint_from_edge(size_t ei) const { return DartConstraints[Edges[ei].Constraint]; }; std::shared_ptr select_contour(Point p); std::shared_ptr contour(size_t i) const { return Contours[i]; }; std::shared_ptr curve_from_edge(size_t ei) const { return DartConstraints[Edges[ei].Constraint].curve(); }; std::shared_ptr curve(size_t i) const {return BoundaryConstraints[i]->curve(); }; std::shared_ptr boundary_constraint(std::shared_ptr const &c) const; std::shared_ptr boundary_constraint(size_t i) const {return BoundaryConstraints[i]; }; std::shared_ptr add_mapped_boundary_pair(ContinuousBoundaryPair const &cbp) { /* * Constructs a MappedBoundaryPair mesh constraint from a ContinuousBoundaryPair * * //TODO: Present implementation strategy allows cyclical mapping * //TODO: which should converge for boundaries restricted to power-of-2 parametric discretizations, but is untested */ MappedBoundaries.push_back(std::make_shared(*this, cbp)); return MappedBoundaries.back(); } template std::vector> add_mapped_boundary_pair(std::vector const &cbp_vec) { /* * Constructs a MappedBoundaryPair for each ContinuousBoundaryPair in cbp_vec */ for(T const &cbp : cbp_vec) { MappedBoundaries.push_back(std::make_shared(*this, cbp)); } return std::vector>(MappedBoundaries.end() - cbp_vec.size(), MappedBoundaries.end()); } DartConstraint const dart_constraint(size_t i) const { return DartConstraints[i]; }; DartTriangle const &triangle(size_t i) const { return Triangles[i]; }; Point circumcenter(size_t ei) const; Point const &base(Edge const &e) const { return Points[e.Node]; }; Point const &base(size_t i) const { return Points[node(i)]; }; Point const &point(size_t i) const { return Points[i]; }; Point const &point(Edge const &e) const { return Points[e.Node]; }; Point const &tip(Edge const &e) const { return Points[next(e).Node]; }; Point const &tip(size_t i) const { return Points[node(next(i))]; }; Point midpoint(size_t e) { return midpoint(Edges[e]); }; Point midpoint(Edge const e) const { Point const &p0 = base(e); Point const &p1 = tip(e); return Point((p0.X + p1.X) / 2.0, (p0.Y + p1.Y) / 2.0); } Edge const &edge(size_t i) const { return Edges[i]; }; Edge const &next(Edge const &e) const { return Edges[e.Next]; }; Edge const &prev(Edge const &e) const { return Edges[e.Prev]; }; Edge const &twin(Edge const &e) const { return Edges[e.Twin]; }; Edge const &edge_from_triangle_index(size_t i) const { return Edges[Triangles[i].Edge]; }; // TODO: Rename to triangle_dart Edge const &edge_from_triangle(DartTriangle const &dt) const { return Edges[dt.Edge]; }; std::array nodes_from_triangle(DartTriangle const & dt) const { Edge const &e = edge_from_triangle(dt); std::array n; n[0] = e.Node; n[1] = next(e).Node; n[2] = prev(e).Node; return n; }; std::vector boundary_nodes(size_t i) const { return BoundaryConstraints[i]->nodes(*this); } LocateTriangleResult locate_triangle(Point const p, size_t &ei) const; LocateTriangleResult locate_triangle(Point const p) const { size_t ei = Edges.size() - 1; return locate_triangle(p, ei); }; AddToQueueResult add_circumcenter_to_queue(size_t dart); AddToQueueResult add_point_to_queue(Point const p, size_t dart); AddToQueueResult add_point_to_queue(Point const p) { return add_point_to_queue(p, Edges.size() - 1); }; protected: std::shared_ptr Boundary; std::vector> Contours; std::vector> BoundaryConstraints; std::vector> MappedBoundaries; std::vector Points; std::vector Edges; std::vector DartConstraints; std::vector Triangles; std::vector> Queue; private: bool find_attached(Point const p, size_t &ei); bool recursive_swap(size_t ei); bool swap(size_t ei); size_t new_edges(size_t num_new) { for (size_t i = 0; i != num_new; ++i) { Edges.emplace_back(Edges.size()); Edges.back().Constraint = 0; } return Edges.size(); } DartConstraint &new_dart_constraint(DartConstraint &dc) { return new_dart_constraint(dc.S0, dc.S1, dc.boundary_constraint()); } DartConstraint &new_dart_constraint(double s0, double s1, std::shared_ptr bc) { bc->add_dart_constraint(DartConstraints.size()); DartConstraints.emplace_back(s0, s1, bc, DartConstraints.size()); return DartConstraints.back(); } Edge &new_edge(size_t p, size_t c, bool dir) { Edges.emplace_back(p, Edges.size(), c, dir); return Edges.back(); } void add_dart_to_queue(size_t d); void add_encroached_edges_to_queue_and_insert(); void create_boundary_polygon(); void element_quality(std::vector &radii, std::vector &quality); void enforce_boundary_mappings(); void get_triangles(); void insert_from_queue(); void insert_internal_boundaries(); void make_edges_optimal(); void mark_triangles(); void refine_once(std::vector index, std::vector circumradius, std::vector quality); void sort_constraints(); void sort_constraints_by_length(); void sort_permutation_ascending(std::vector &value, std::vector &index) const; void sort_permutation_descending(std::vector &values, std::vector &index) const; void split_edge(size_t ei); void split_encroached_edges(); void triangulate_boundary_polygon(); std::vector get_encroached_edges(Point const p, size_t edge); InsertPointResult insert_point(Point const p, size_t ei); InsertPointResult insert_midpoint(size_t ei); }; #endif //OERSTED_MESH_H