BoundaryCondition.h 5.56 KB
Newer Older
JasonPries's avatar
JasonPries committed
1 2 3 4 5 6 7 8 9
#ifndef OERSTED_BOUNDARYCONDITION_H
#define OERSTED_BOUNDARYCONDITION_H

#include <vector>
#include <set>
#include <algorithm>

#include <Eigen/Sparse>

10 11
#include "Quadrature.hpp"

JasonPries's avatar
JasonPries committed
12 13 14
#include "MatrixGroup.h"
#include "FiniteElementMesh.h"

JasonPries's avatar
JasonPries committed
15 16 17
// TODO: DiscreteBoundaryCondition versus ContinuousBoundaryCondition? Somehow need to apply to top level model and have that be translated to the discrete model
// TODO: Boundary conditions should represent simple containers that know how to perform actions on matrices, parsing of boundaries to extract nodes should take place at a higher level

JasonPries's avatar
JasonPries committed
18 19
class BoundaryCondition {
public:
JasonPries's avatar
JasonPries committed
20
    virtual ~BoundaryCondition() {};
JasonPries's avatar
JasonPries committed
21 22

    virtual void apply(std::vector<Eigen::Triplet<double>> &triplets) const = 0;
JasonPries's avatar
JasonPries committed
23 24

    virtual void reduce(std::set<size_t, std::less<size_t>> &index) const = 0;
JasonPries's avatar
JasonPries committed
25 26 27
};

template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
28 29
class ZeroDirichlet : public BoundaryCondition {
};
JasonPries's avatar
JasonPries committed
30 31

template<size_t ElementOrder, size_t QuadratureOrder>
JasonPries's avatar
JasonPries committed
32
class ZeroDirichlet<2, ElementOrder, QuadratureOrder> : public BoundaryCondition {
JasonPries's avatar
JasonPries committed
33
public:
JasonPries's avatar
JasonPries committed
34
    ZeroDirichlet(std::vector<std::shared_ptr<Boundary<2>>> boundaries) : Boundaries{boundaries} {};
JasonPries's avatar
JasonPries committed
35

36 37 38 39 40 41 42
    ZeroDirichlet(std::vector<std::shared_ptr<Curve const>> const &boundaries, FiniteElementMesh<2, ElementOrder> const &fem) {
        for (auto const &b : boundaries) {
            auto fb = fem.boundary(b);
            if (fb) { Boundaries.push_back(fb); }
        }
    }

JasonPries's avatar
JasonPries committed
43 44
    void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {};

JasonPries's avatar
JasonPries committed
45 46 47
    void reduce(std::set<size_t, std::less<size_t>> &index) const override {
        for (auto const &b : Boundaries) {
            for (auto const &i : b->nodes()) {
JasonPries's avatar
JasonPries committed
48 49
                index.insert(i);
            }
JasonPries's avatar
JasonPries committed
50 51 52 53
        }
    };

protected:
JasonPries's avatar
JasonPries committed
54
    std::vector<std::shared_ptr<Boundary<2>>> Boundaries;
JasonPries's avatar
JasonPries committed
55
};
JasonPries's avatar
JasonPries committed
56

JasonPries's avatar
JasonPries committed
57
// TODO: Unify interface between PeriodicBoundary and SlidingBoundary (inherit a common interface/data members)
JasonPries's avatar
JasonPries committed
58 59

template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
60 61
class PeriodicBoundaryCondition : public BoundaryCondition {
};
JasonPries's avatar
JasonPries committed
62 63

template<size_t ElementOrder, size_t QuadratureOrder>
JasonPries's avatar
JasonPries committed
64
class PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder> : public BoundaryCondition {
JasonPries's avatar
JasonPries committed
65
public:
66
    PeriodicBoundaryCondition(std::vector<VariableMap> map, bool antiperiodic) : Map{map}, Antiperiodic{antiperiodic} {};
JasonPries's avatar
JasonPries committed
67

68 69
    PeriodicBoundaryCondition(std::vector<PeriodicBoundaryPair> const &periodic_boundary_pairs, FiniteElementMesh<2, ElementOrder> const &fem, bool antiperiodic) : Antiperiodic{antiperiodic} {
        double_t sgn = sign();
70

71
        for (auto const &pbp : periodic_boundary_pairs) {
72 73 74 75 76
            auto const &n0 = fem.boundary(pbp.curve0())->nodes();
            auto const &n1 = fem.boundary(pbp.curve1())->nodes();

            if (pbp.orientation()) {
                for (size_t i = 0; i != n0.size(); ++i) {
77
                    Map.emplace_back(n0[i], n1[i], sgn);
78 79 80 81
                }
            } else {
                size_t j = n1.size();
                for (size_t i = 0; i != n0.size(); ++i) {
82
                    Map.emplace_back(n0[i], n1[--j], sgn);
83 84
                }
            }
85 86 87 88 89 90
        }

        // Remove duplicates
        std::sort(Map.begin(), Map.end());
        auto last = std::unique(Map.begin(), Map.end());
        Map.erase(last, Map.end());
91

92 93 94 95 96 97
        if (!antiperiodic) { // Remove center node if it exists
            auto eq = [](VariableMap const &x) { return (x.first() == x.second()); };
            auto center = std::find_if(Map.begin(), Map.end(), eq);
            if (center != Map.end()) { Map.erase(center); };
        }
        // TODO: needs testing
98 99
    }

JasonPries's avatar
JasonPries committed
100
    void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {
101
        //double_t sgn = sign();
JasonPries's avatar
JasonPries committed
102

103
        for (VariableMap const &v : Map) {
104
            triplets.push_back(Eigen::Triplet<double>(v.first(), v.second(), v.value()));
JasonPries's avatar
JasonPries committed
105
        }
JasonPries's avatar
JasonPries committed
106 107
    }

JasonPries's avatar
JasonPries committed
108
    void reduce(std::set<size_t, std::less<size_t>> &index) const override {
109 110
        for (VariableMap const &v : Map) {
            index.insert(v.second());
JasonPries's avatar
JasonPries committed
111 112 113
        }
    }

114 115
    double_t sign() const { return (Antiperiodic ? -1.0 : 1.0); };

JasonPries's avatar
JasonPries committed
116
protected:
117
    std::vector<VariableMap> Map;
JasonPries's avatar
JasonPries committed
118 119 120 121

    bool Antiperiodic;
};

JasonPries's avatar
JasonPries committed
122
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
123 124
class SlidingInterface : public BoundaryCondition {
};
JasonPries's avatar
JasonPries committed
125

JasonPries's avatar
JasonPries committed
126
template<size_t ElementOrder, size_t QuadratureOrder>
127
class SlidingInterface<2, ElementOrder, QuadratureOrder> : public BoundaryCondition {
JasonPries's avatar
JasonPries committed
128 129 130 131
public:
    SlidingInterface(std::vector<size_t> first_vars, std::vector<size_t> second_vars, std::vector<double> vals, bool antiperiodic = false)
            : First{first_vars}, Second{second_vars}, Value{vals}, Antiperiodic{antiperiodic} {};

132
    SlidingInterface &operator++() {
JasonPries's avatar
JasonPries committed
133 134 135 136 137 138 139
        std::rotate(Second.begin(), Second.begin() + 1, Second.end());
        if (Antiperiodic) {
            std::rotate(Value.begin(), Value.begin() + 1, Value.end());
            *(Value.end()) = -*(Value.end());
        }
    }

JasonPries's avatar
JasonPries committed
140
    // TODO: SlidingInterface& operator--() {...};
JasonPries's avatar
JasonPries committed
141 142

    void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {
143 144
        for (size_t i = 0; i != Value.size(); ++i) {
            triplets.push_back(Eigen::Triplet<double>(First[i], Second[i], Value[i]));
JasonPries's avatar
JasonPries committed
145 146 147 148
        }
    }

    void reduce(std::set<size_t, std::less<size_t>> &index) const override {
149
        for (size_t i = 0; i != Second.size(); ++i) {
JasonPries's avatar
JasonPries committed
150 151 152 153 154 155 156 157 158 159 160
            index.insert(Second[i]);
        }
    }

protected:
    std::vector<size_t> First;
    std::vector<size_t> Second;
    std::vector<double> Value;

    bool Antiperiodic;
};
161

JasonPries's avatar
JasonPries committed
162
#endif //OERSTED_BOUNDARYCONDITION_H