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

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

#include <Eigen/Sparse>

#include "MatrixGroup.h"
#include "QuadratureRule.h"
#include "FiniteElementMesh.h"

JasonPries's avatar
JasonPries committed
14
15
16
// 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
17
18
class BoundaryCondition {
public:
JasonPries's avatar
JasonPries committed
19
    virtual ~BoundaryCondition() {};
JasonPries's avatar
JasonPries committed
20
21

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

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

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

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

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

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

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

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

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

template<size_t ElementOrder, size_t QuadratureOrder>
JasonPries's avatar
JasonPries committed
54
class PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder> : public BoundaryCondition {
JasonPries's avatar
JasonPries committed
55
public:
JasonPries's avatar
JasonPries committed
56
    PeriodicBoundaryCondition(std::vector<BoundaryMap> map, bool antiperiodic) : Map{map}, Antiperiodic{antiperiodic} {};
JasonPries's avatar
JasonPries committed
57
58
59
60

    void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {
        double sign{1.0 - 2.0 * Antiperiodic};

JasonPries's avatar
JasonPries committed
61
62
63
64
65
        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
66
67
    }

JasonPries's avatar
JasonPries committed
68
69
70
71
    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
72
            }
JasonPries's avatar
JasonPries committed
73
74
75
76
        }
    }

protected:
JasonPries's avatar
JasonPries committed
77
    std::vector<BoundaryMap> Map;
JasonPries's avatar
JasonPries committed
78
79
80
81

    bool Antiperiodic;
};

JasonPries's avatar
JasonPries committed
82
83
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
class SlidingInterface : public BoundaryCondition {};
JasonPries's avatar
JasonPries committed
84

JasonPries's avatar
JasonPries committed
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
112
113
114
115
116
117
118
119
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());
        }
    }

    // TODO: SlidingInterface& opertor--() {...};

    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
120
#endif //OERSTED_BOUNDARYCONDITION_H