BoundaryCondition.h 5.35 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>
JasonPries's avatar
JasonPries committed
28
class ZeroDirichlet : public BoundaryCondition {};
JasonPries's avatar
JasonPries committed
29
30

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

    void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {};

JasonPries's avatar
JasonPries committed
37
38
39
    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
40
41
                index.insert(i);
            }
JasonPries's avatar
JasonPries committed
42
43
44
45
        }
    };

protected:
JasonPries's avatar
JasonPries committed
46
    std::vector<std::shared_ptr<Boundary<2>>> Boundaries;
JasonPries's avatar
JasonPries committed
47
};
JasonPries's avatar
JasonPries committed
48

JasonPries's avatar
JasonPries committed
49
// TODO: Unify interface between PeriodicBoundary and SlidingBoundary (inherit a common interface/data members)
JasonPries's avatar
JasonPries committed
50
51

template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
JasonPries's avatar
JasonPries committed
52
class PeriodicBoundaryCondition : public BoundaryCondition {};
JasonPries's avatar
JasonPries committed
53
54

template<size_t ElementOrder, size_t QuadratureOrder>
JasonPries's avatar
JasonPries committed
55
class PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder> : public BoundaryCondition {
JasonPries's avatar
JasonPries committed
56
public:
JasonPries's avatar
JasonPries committed
57
    PeriodicBoundaryCondition(std::vector<BoundaryMap> map, bool antiperiodic) : Map{map}, Antiperiodic{antiperiodic} {};
JasonPries's avatar
JasonPries committed
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
    PeriodicBoundaryCondition(std::vector<PeriodicBoundaryPair> const &periodic_boundary_pairs, FiniteElementMesh<2,ElementOrder> const &fem, bool antiperiodic) : Antiperiodic{antiperiodic} {
        Map.reserve(periodic_boundary_pairs.size());

        double_t sign = (antiperiodic ? -1.0 : 1.0);
        for(auto const &pbp : periodic_boundary_pairs) {
            // TODO: Change Map type to std::set<VariableMap>, add comparison operator to VariableMap, and make sure node indicies are lexigraphically ordered in VariableMap
            auto const &n0 = fem.boundary(pbp.curve0())->nodes();
            auto const &n1 = fem.boundary(pbp.curve1())->nodes();
            std::vector<VariableMap> vmap;
            vmap.reserve(n0.size());

            if (pbp.orientation()) {
                for (size_t i = 0; i != n0.size(); ++i) {
                    vmap.emplace_back(n0[i], n1[i], sign);
                }
            } else {
                size_t j = n1.size();
                for (size_t i = 0; i != n0.size(); ++i) {
                    vmap.emplace_back(n0[i], n1[--j], sign);
                }
            }

            Map.emplace_back(vmap);// TODO: needs testing, allows duplicate nodes which may not be handle correctly by constraint algorithm
         }
    }

JasonPries's avatar
JasonPries committed
85
86
87
    void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {
        double sign{1.0 - 2.0 * Antiperiodic};

JasonPries's avatar
JasonPries committed
88
89
90
91
92
        for (BoundaryMap const &b : Map) {
            for (VariableMap const &v : b.map()) {
                triplets.push_back(Eigen::Triplet<double>(v.first(), v.second(), sign * v.value()));
            }
        }
JasonPries's avatar
JasonPries committed
93
94
    }

JasonPries's avatar
JasonPries committed
95
96
97
98
    void reduce(std::set<size_t, std::less<size_t>> &index) const override {
        for (BoundaryMap const &b : Map) {
            for (VariableMap const &v : b.map()) {
                    index.insert(v.second());
JasonPries's avatar
JasonPries committed
99
            }
JasonPries's avatar
JasonPries committed
100
101
102
103
        }
    }

protected:
JasonPries's avatar
JasonPries committed
104
    std::vector<BoundaryMap> Map;
JasonPries's avatar
JasonPries committed
105
106
107
108

    bool Antiperiodic;
};

JasonPries's avatar
JasonPries committed
109
110
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
class SlidingInterface : public BoundaryCondition {};
JasonPries's avatar
JasonPries committed
111

JasonPries's avatar
JasonPries committed
112
113
114
115
116
117
118
119
120
121
122
123
124
125
template<size_t ElementOrder, size_t QuadratureOrder>
class SlidingInterface<2,ElementOrder,QuadratureOrder> : public BoundaryCondition {
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} {};

    SlidingInterface& operator++() {
        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
126
    // TODO: SlidingInterface& operator--() {...};
JasonPries's avatar
JasonPries committed
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146

    void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {
        for(size_t i = 0; i!=Value.size();++i){
            triplets.push_back(Eigen::Triplet<double>(First[i],Second[i],Value[i]));
        }
    }

    void reduce(std::set<size_t, std::less<size_t>> &index) const override {
        for (size_t i = 0; i!=Second.size();++i){
            index.insert(Second[i]);
        }
    }

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

    bool Antiperiodic;
};
JasonPries's avatar
JasonPries committed
147
#endif //OERSTED_BOUNDARYCONDITION_H