BoundaryCondition.h

#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"

// 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

class BoundaryCondition {
public:
    virtual ~BoundaryCondition() {};

    virtual void apply(std::vector<Eigen::Triplet<double>> &triplets) const = 0;

    virtual void reduce(std::set<size_t, std::less<size_t>> &index) const = 0;
};

template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
class ZeroDirichlet : public BoundaryCondition {};

template<size_t ElementOrder, size_t QuadratureOrder>
class ZeroDirichlet<2, ElementOrder, QuadratureOrder> : public BoundaryCondition {
public:
    ZeroDirichlet(std::vector<std::shared_ptr<Boundary<2>>> boundaries) : Boundaries{boundaries} {};

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

    void reduce(std::set<size_t, std::less<size_t>> &index) const override {
        for (auto const &b : Boundaries) {
            for (auto const &i : b->nodes()) {
                index.insert(i);
            }
        }
    };

protected:
    std::vector<std::shared_ptr<Boundary<2>>> Boundaries;
};

// TODO: Unify interface between PeriodicBoundary and SlidingBoundary (inherit a common interface/data members)

template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
class PeriodicBoundaryCondition : public BoundaryCondition {};

template<size_t ElementOrder, size_t QuadratureOrder>
class PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder> : public BoundaryCondition {
public:
    PeriodicBoundaryCondition(std::vector<BoundaryMap> map, bool antiperiodic) : Map{map}, Antiperiodic{antiperiodic} {};

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

        for (BoundaryMap const &b : Map) {
            for (VariableMap const &v : b.map()) {
                triplets.push_back(Eigen::Triplet<double>(v.first(), v.second(), sign * v.value()));
            }
        }
    }

    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());
            }
        }
    }

protected:
    std::vector<BoundaryMap> Map;

    bool Antiperiodic;
};

template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
class SlidingInterface : public BoundaryCondition {};

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& operator--() {...};

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