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