Mesh.h 9.27 KB
Newer Older
1 2 3 4 5 6
#ifndef OERSTED_MESH_H
#define OERSTED_MESH_H

#include <algorithm>
#include <fstream>
#include <numeric>
7
#include <set>
8

9 10 11 12
#include "Eigen"
#include "Eigen/Sparse"
#include "Eigen/SparseLU"

13 14
#include "Sketch.hpp"

15
#include "BoundaryConstraint.h"
16
#include "MappedBoundaryPair.h"
17 18
#include "Edge.h"
#include "InsertionQueuer.h"
19
#include "Point.h"
20
#include "DartTriangle.h"
21

22 23 24 25 26 27 28 29
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
};

30
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)
31
public:
JasonPries's avatar
JasonPries committed
32 33
    Mesh() {};

34 35
    Mesh(Sketch &s);

36
    friend class BoundaryConstraint;
37
    friend class MappedBoundaryPair;
38 39 40
    friend class CircumcenterQueuer;
    friend class MidpointQueuer;

41 42 43 44 45 46
    double MinimumElementQuality = 0.0;
    double MinimumElementSize = 0.0;
    double MaximumElementSize = DBL_MAX;

    bool are_intersecting(size_t ei, size_t ej) const;

47 48
    bool edges_are_optimal() const;

49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
    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; };

89
    size_t size_dart_constraints() const {return DartConstraints.size(); };
90 91 92

    size_t size_edges() const { return Edges.size(); };

93 94
    size_t size_points() const { return Points.size(); };

95 96
    size_t size_triangles() const { return Triangles.size(); };

97 98
    size_t size_contours() const { return Contours.size(); };

99 100
    size_t size_boundary_constraints() const { return BoundaryConstraints.size(); };

101 102
    size_t twin(size_t ei) const { return Edges[ei].Twin; };

103 104
    void add_to_queue(std::unique_ptr<InsertionQueuer> ic) { Queue.push_back(std::move(ic)); };

105 106 107 108
    void create();

    void save_as(std::string path, std::string file_name) const;

109 110
    void smooth();

111
    // TODO: rename these methods something like, dart_constraint_from_edge_index, curve_from_edge_index
112 113
    DartConstraint const dart_constraint_from_edge(size_t ei) const { return DartConstraints[Edges[ei].Constraint]; };

114 115 116 117
    std::shared_ptr<Contour const> select_contour(Point p);

    std::shared_ptr<Contour const> contour(size_t i) const { return Contours[i]; };

118
    std::shared_ptr<Curve const> curve_from_edge(size_t ei) const { return DartConstraints[Edges[ei].Constraint].curve(); };
119

120
    std::shared_ptr<Curve const> curve(size_t i) const {return BoundaryConstraints[i]->curve(); };
121

122 123
    std::shared_ptr<BoundaryConstraint> boundary_constraint(std::shared_ptr<Curve const> const &c) const;

124 125
    std::shared_ptr<BoundaryConstraint> boundary_constraint(size_t i) const {return BoundaryConstraints[i]; };

126 127 128 129 130 131 132 133 134 135 136 137
    std::shared_ptr<MappedBoundaryPair> 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<MappedBoundaryPair>(*this, cbp));
        return MappedBoundaries.back();
    }

138 139 140 141 142 143 144 145 146 147 148 149 150
    template<class T>
    std::vector<std::shared_ptr<MappedBoundaryPair>> add_mapped_boundary_pair(std::vector<T> const &cbp_vec) {
        /*
         * Constructs a MappedBoundaryPair for each ContinuousBoundaryPair in cbp_vec
         */

        for(T const &cbp : cbp_vec) {
            MappedBoundaries.push_back(std::make_shared<MappedBoundaryPair>(*this, cbp));
        }

        return std::vector<std::shared_ptr<MappedBoundaryPair>>(MappedBoundaries.end() - cbp_vec.size(), MappedBoundaries.end());
    }

151
    DartConstraint const dart_constraint(size_t i) const { return DartConstraints[i]; };
152

153 154 155 156
    DartTriangle const &triangle(size_t i) const { return Triangles[i]; };

    Point circumcenter(size_t ei) const;

157
    Point const &base(Edge const &e) const { return Points[e.Node]; };
158

159
    Point const &base(size_t i) const { return Points[node(i)]; };
160

161
    Point const &point(size_t i) const { return Points[i]; };
162

163
    Point const &point(Edge const &e) const { return Points[e.Node]; };
164

165
    Point const &tip(Edge const &e) const { return Points[next(e).Node]; };
166

167
    Point const &tip(size_t i) const { return Points[node(next(i))]; };
168

169
    Point midpoint(size_t e) { return midpoint(Edges[e]); };
170

171 172 173
    Point midpoint(Edge const e) const {
        Point const &p0 = base(e);
        Point const &p1 = tip(e);
174

175 176 177 178 179 180 181 182
        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]; };
183

184
    Edge const &twin(Edge const &e) const { return Edges[e.Twin]; };
185

186 187
    Edge const &edge_from_triangle_index(size_t i) const { return Edges[Triangles[i].Edge]; }; // TODO: Rename to triangle_dart

188
    Edge const &edge_from_triangle(DartTriangle const &dt) const { return Edges[dt.Edge]; };
189

190
    std::array<size_t,3> nodes_from_triangle(DartTriangle const & dt) const {
191 192 193 194 195 196 197 198 199 200
        Edge const &e = edge_from_triangle(dt);

        std::array<size_t,3> n;

        n[0] = e.Node;
        n[1] = next(e).Node;
        n[2] = prev(e).Node;

        return n;
    };
201

202 203
    std::vector<size_t> boundary_nodes(size_t i) const { return BoundaryConstraints[i]->nodes(*this); }

204 205 206 207 208 209 210
    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);
    };

211 212 213 214 215
    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); };
216 217 218 219

protected:
    std::shared_ptr<Contour const> Boundary;
    std::vector<std::shared_ptr<Contour const>> Contours;
220
    std::vector<std::shared_ptr<BoundaryConstraint>> BoundaryConstraints;
221
    std::vector<std::shared_ptr<MappedBoundaryPair>> MappedBoundaries;
222 223 224

    std::vector<Point> Points;
    std::vector<Edge> Edges;
225
    std::vector<DartConstraint> DartConstraints;
226
    std::vector<DartTriangle> Triangles;
227

228
    std::vector<std::unique_ptr<InsertionQueuer>> Queue;
229

230 231 232 233 234 235 236 237 238
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) {
239
            Edges.emplace_back(Edges.size());
240 241 242 243 244
            Edges.back().Constraint = 0;
        }
        return Edges.size();
    }

245 246 247 248 249
    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<BoundaryConstraint> bc) {
250
        bc->add_dart_constraint(DartConstraints.size());
251 252 253 254
        DartConstraints.emplace_back(s0, s1, bc, DartConstraints.size());
        return DartConstraints.back();
    }

255
    Edge &new_edge(size_t p, size_t c, bool dir) {
256
        Edges.emplace_back(p, Edges.size(), c, dir);
257 258 259
        return Edges.back();
    }

260 261 262 263
    void add_dart_to_queue(size_t d);

    void add_encroached_edges_to_queue_and_insert();

264 265 266 267
    void create_boundary_polygon();

    void element_quality(std::vector<double> &radii, std::vector<double> &quality);

268 269
    void enforce_boundary_mappings();

270 271
    void get_triangles();

272 273
    void insert_from_queue();

274 275
    void insert_internal_boundaries();

276 277
    void make_edges_optimal();

278 279 280 281
    void mark_triangles();

    void refine_once(std::vector<size_t> index, std::vector<double> circumradius, std::vector<double> quality);

282 283
    void sort_constraints();

284 285
    void sort_constraints_by_length();

286 287 288 289 290 291 292 293 294 295
    void sort_permutation_ascending(std::vector<double> &value, std::vector<size_t> &index) const;

    void sort_permutation_descending(std::vector<double> &values, std::vector<size_t> &index) const;

    void split_edge(size_t ei);

    void split_encroached_edges();

    void triangulate_boundary_polygon();

296
    std::vector<size_t> get_encroached_edges(Point const p, size_t edge);
297 298 299 300 301 302 303

    InsertPointResult insert_point(Point const p, size_t ei);

    InsertPointResult insert_midpoint(size_t ei);
};

#endif //OERSTED_MESH_H