BoundaryConstraint.cpp 4.93 KB
Newer Older
1 2
#include "BoundaryConstraint.h"
#include "DartConstraint.h"
3
#include "Mesh.h"
4

5 6 7 8
bool BoundaryConstraint::is_uniform(Mesh const &m) const {
    DartConstraint const &dc = m.dart_constraint(DartConstraints[0]);
    double_t delta_0{dc.delta()};

9
    for (size_t i : DartConstraints) {
10 11 12
        double delta_i{m.dart_constraint(i).delta()};
        if (std::abs(delta_0 - delta_i) > FLT_EPSILON) {
            return false;
13 14
        }
    }
15 16

    return true;
17 18
}

19 20
double BoundaryConstraint::smallest_parametric_edge_length(Mesh &m) const {
    double_t delta_min{DBL_MAX};
21

22 23 24 25
    for (size_t i : DartConstraints) {
        DartConstraint const &dc = m.dart_constraint(i);
        delta_min = std::min(delta_min, std::abs(dc.S1 - dc.S0));
    }
26

27 28
    return delta_min;
}
29

30
double BoundaryConstraint::queue_uniform(Mesh &m, double delta_min) const {
31
    std::cout << "//TODO: Need to make sure that all edges encroached by the edge mid-point are enqueued" << std::endl;
32
    std::cout << "//TODO: Constrained dart enqueuing needs to be revisited, and need to make sure enqueued edges are unique" << std::endl;
33
    double delta_max{delta_min};
34

35
    for (size_t i : DartConstraints) {
36 37 38
        DartConstraint const &dc = m.dart_constraint(i);
        double_t delta = std::abs(dc.S1 - dc.S0);
        if (delta > delta_min) {
39 40
            //m.add_to_queue(std::make_unique<MidpointQueuer>(dc.dart()));
            m.add_dart_to_queue(dc.dart());
41
            delta_max = std::max(delta_max, delta);
42
        }
43
    }
44

45 46 47 48 49 50 51
    return delta_max;
}

void BoundaryConstraint::add_to_queue(Mesh &m, size_t dart) const {
    if (UniformDiscretization) {
        queue_uniform(m, smallest_parametric_edge_length(m) / 2.0);
    } else {
52 53
        //m.add_dart_to_queue(dart);
        m.add_to_queue(std::make_unique<MidpointQueuer>(dart));
54 55 56
    }
}

57 58 59 60 61 62
void BoundaryConstraint::insert_parameters(Mesh &mesh, std::vector<double_t> s) {
    while (s.size() > 0) {
        queue_parameters(mesh, s);
        mesh.insert_from_queue();
    }
}
63

64
void BoundaryConstraint::make_uniform(Mesh &m) { make_uniform(m, smallest_parametric_edge_length(m)); };
65

66 67
void BoundaryConstraint::make_uniform(Mesh &m, double delta_min) {
    delta_min *= (1 + FLT_EPSILON);
68 69
    double_t delta_max{DBL_MAX};
    while (delta_max > delta_min) {
70 71 72 73
        delta_max = queue_uniform(m, delta_min);
        m.insert_from_queue();
    }
}
74

75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
void BoundaryConstraint::queue_parameters(Mesh &mesh, std::vector<double_t> &s) {
    /*
        Enqueues all darts having midpoints coincident with the parameter vector s
        The parameter vector s is consumed, with s[i] being erased if an associated dart is enqueued

        TODO: This method assumes all values in s are powers of 2 and should check that assumption
     */

    sort_dart_constraints(mesh);
    std::sort(s.begin(), s.end());

    size_t j{0};
    for (size_t i : DartConstraints) {
        auto const &dc = mesh.dart_constraint(i);

        while (j != s.size() && s[j] < dc.S1 + FLT_EPSILON) {
            if (std::abs(dc.S0 + dc.S1 - 2.0 * s[j]) < FLT_EPSILON) { // if s[j] is the dart midpoint, split the dart in half
                mesh.add_dart_to_queue(dc.dart());

                s.erase(s.begin() + j);
                break;
            }

            if (std::abs(dc.S0 - s[j]) < FLT_EPSILON || std::abs(dc.S1 - s[j]) < FLT_EPSILON) { // if s[j] is either dart endpoint, remove s from the queue
                s.erase(s.begin() + j);
                break;
            }

            ++j;
        }

        if (s.size() == 0) {
            break;
        }
    }
}

112 113 114 115 116 117
void BoundaryConstraint::refine(Mesh &m, double_t tol) {
    double_t err_max{DBL_MAX};

    while (err_max > tol) {
        err_max = 0.0;
        for (size_t i : DartConstraints) {
118
            DartConstraint const &dc = m.dart_constraint(i);
119 120 121 122 123 124

            double_t len_c = ConstraintCurve->length(dc.S0, dc.S1);
            double_t len_e = m.edge(dc.dart()).length(m);
            double_t err_i = std::abs(len_c - len_e) / len_c;

            if (err_i > tol) {
125 126
                m.add_to_queue(std::make_unique<MidpointQueuer>(dc.dart()));

127
                err_max = std::max(err_max, err_i);
128 129 130
            }
        }

131 132
        m.insert_from_queue();
    }
133 134 135
}

void BoundaryConstraint::sort_dart_constraints(Mesh const &m) {
136 137
    auto comp = [&](size_t i, size_t j) { return m.dart_constraint(i).S0 < m.dart_constraint(j).S0; };
    std::sort(DartConstraints.begin(), DartConstraints.end(), comp);
138 139 140 141 142 143 144 145 146
}

std::vector<size_t> BoundaryConstraint::nodes(Mesh const &m) {
    sort_dart_constraints(m);
    std::vector<size_t> n;
    n.reserve(DartConstraints.size());
    for (size_t i : DartConstraints) {
        n.push_back(m.dart_constraint(i).base(m));
    }
147
    n.push_back(m.dart_constraint(DartConstraints.back()).tip(m));
148 149

    return n;
150 151 152 153 154 155 156 157 158 159
}

std::vector<double_t> BoundaryConstraint::parameters(Mesh const &m) {
    std::vector<double_t> s;
    for (size_t d : DartConstraints) {
        s.emplace_back(m.dart_constraint(d).S0);
    }
    s.emplace_back(1.0);

    return s;
160
}