BoundaryCondition.h 4.1 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

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

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

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

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

    bool Antiperiodic;
};

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

JasonPries's avatar
JasonPries committed
86
87
88
89
90
91
92
93
94
95
96
97
98
99
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
100
    // TODO: SlidingInterface& operator--() {...};
JasonPries's avatar
JasonPries committed
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120

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