Physics.h 10.3 KB
Newer Older
1 2 3 4
#ifndef OERSTED_PHYSICS_H
#define OERSTED_PHYSICS_H

#include <memory>
JasonPries's avatar
JasonPries committed
5
#include <set>
6

JasonPries's avatar
JasonPries committed
7
#include "Eigen/IterativeLinearSolvers"
JasonPries's avatar
JasonPries committed
8
#include "Eigen/SparseCholesky"
JasonPries's avatar
JasonPries committed
9

JasonPries's avatar
JasonPries committed
10
#include "BoundaryCondition.h"
11
#include "FiniteElementMesh.h"
JasonPries's avatar
JasonPries committed
12
#include "Forcing.h"
JasonPries's avatar
JasonPries committed
13 14 15
#include "Node.h"
#include "PhysicsConstants.h"
#include "Triangle.h"
16

JasonPries's avatar
JasonPries committed
17 18 19 20
// TODO: These classes should be generic (e.g. CurlCurl, ScalarLaplacian, VectorLaplacian)
// TODO:    Higher level classes should implement interfaces with appropriate physics grammar
// TODO:    Implement material property dependencies using pointers to virtual member functions

21
enum class FieldVariable {
22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
    MagneticVectorPotential = 1,
    A = 1,
    MagneticFluxDensity = 2,
    B = 2,
    ElectricFluxDensity = 4,
    D = 4,
    ElectricFieldIntensity = 5,
    E = 5,
    MagneticFieldIntensity = 8,
    H = 8,
    ElectricCurrentDensity = 10,
    J = 10,
    ElectricScalarPotential = 26 + 21,
    Phi = 26 + 21,
    MagneticScalarPotential = 26 + 23,
    Psi = 26 + 23
};

JasonPries's avatar
JasonPries committed
40 41 42 43 44 45 46 47 48
class PhysicsInterface {
public:
    double time() { return Time; };

    void time(double t) { Time = t; };

protected:
    double Time;

JasonPries's avatar
JasonPries committed
49 50
    std::vector<std::shared_ptr<Forcing>> ForcingCondtions;
    std::vector<std::shared_ptr<BoundaryCondition>> BoundaryConditions;
JasonPries's avatar
JasonPries committed
51 52
};

JasonPries's avatar
JasonPries committed
53
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder>
JasonPries's avatar
JasonPries committed
54
class PhysicsData {
JasonPries's avatar
JasonPries committed
55
public:
56
    static constexpr size_t QuadraturePoints = TriangleQuadrature<QuadratureOrder>::size;
JasonPries's avatar
JasonPries committed
57

JasonPries's avatar
JasonPries committed
58
    /* TODO: template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder> class PhysicsData : public PhysicsInterface {...};
JasonPries's avatar
JasonPries committed
59 60 61 62 63 64 65
    virtual Eigen::VectorXd init_residual() = 0;
    virtual Eigen::SparseMatrix<double> init_jacobian() = 0;
    virtual void residual(Eigen::SparseMatrix<double> &J, Eigen::VectorXd &r) = 0;
    virtual void build_matrices() = 0;
    virtual void apply_conditions() = 0;
     */

JasonPries's avatar
JasonPries committed
66
    PhysicsData(FiniteElementMesh<2, ElementOrder> &d) : Domain{d},
JasonPries's avatar
JasonPries committed
67 68 69
                                                         ElementWeights{d.triangles().size()},
                                                         Derivatives{d.nodes().size(), d.triangles().size(), Triangle<ElementOrder>::NumNodes},
                                                         Basis{d.nodes().size(), d.triangles().size(), Triangle<ElementOrder>::NumNodes} {};
JasonPries's avatar
JasonPries committed
70

JasonPries's avatar
JasonPries committed
71
    FiniteElementMesh<Dimension, ElementOrder> const &domain() const { return Domain; };
JasonPries's avatar
JasonPries committed
72 73 74 75 76 77 78 79 80 81 82 83 84

    DiagonalMatrixGroup<QuadraturePoints> const &weights() { return ElementWeights; };

    SparseMatrixGroup<QuadraturePoints> const &basis() { return Basis; };

    DerivativeMatrixGroup<QuadraturePoints> const &derivatives() { return Derivatives; };

protected:
    FiniteElementMesh<Dimension, ElementOrder> &Domain;

    DiagonalMatrixGroup<QuadraturePoints> ElementWeights;
    SparseMatrixGroup<QuadraturePoints> Basis;
    DerivativeMatrixGroup<QuadraturePoints> Derivatives;
85 86
};

87
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder, FieldVariable V>
JasonPries's avatar
JasonPries committed
88
class Magnetostatic : public PhysicsData<Dimension, ElementOrder, QuadratureOrder>, public PhysicsInterface {
JasonPries's avatar
JasonPries committed
89 90 91
};

template<size_t ElementOrder, size_t QuadratureOrder>
92
class Magnetostatic<2, ElementOrder, QuadratureOrder, FieldVariable::MagneticVectorPotential> : public PhysicsData<2, ElementOrder, QuadratureOrder>, public PhysicsInterface {
93
public:
JasonPries's avatar
JasonPries committed
94
    using PhysicsData<2, ElementOrder, QuadratureOrder>::QuadraturePoints;
95

JasonPries's avatar
JasonPries committed
96 97 98 99
    using PhysicsData<2, ElementOrder, QuadratureOrder>::domain;
    using PhysicsData<2, ElementOrder, QuadratureOrder>::derivatives;
    using PhysicsData<2, ElementOrder, QuadratureOrder>::weights;
    using PhysicsData<2, ElementOrder, QuadratureOrder>::basis;
JasonPries's avatar
JasonPries committed
100

JasonPries's avatar
JasonPries committed
101
    Magnetostatic(FiniteElementMesh<2, ElementOrder> &d) : PhysicsData<2, ElementOrder, QuadratureOrder>{d}, Unknowns{d.size_nodes()}, Elements{d.size_elements()} {};
102

JasonPries's avatar
JasonPries committed
103
    Eigen::SparseMatrix<double> init_unknown_matrix() const { return Eigen::SparseMatrix<double>(Unknowns, Unknowns); };
104

JasonPries's avatar
JasonPries committed
105
    Eigen::DiagonalMatrix<double, Eigen::Dynamic> init_element_matrix() const { return Eigen::DiagonalMatrix<double, Eigen::Dynamic>(Elements); }
106

JasonPries's avatar
JasonPries committed
107 108 109
    Eigen::VectorXd init_unknown_vector() const { return Eigen::VectorXd(Unknowns); };

    Eigen::VectorXd init_element_vector() const { return Eigen::VectorXd(Elements); };
110

JasonPries's avatar
JasonPries committed
111 112 113 114
    Eigen::ArrayXd init_unknown_array() const { return Eigen::ArrayXd(Unknowns); };

    Eigen::ArrayXd init_element_array() const { return Eigen::ArrayXd(Elements); };

JasonPries's avatar
JasonPries committed
115
    void linearize(Eigen::SparseMatrix<double> &J, Eigen::VectorXd &v, Eigen::VectorXd &r, Eigen::VectorXd &f,
JasonPries's avatar
JasonPries committed
116
                   Eigen::ArrayXd &Fx, Eigen::ArrayXd &Fy, Eigen::ArrayXd &dFxdx, Eigen::ArrayXd &dFydy, Eigen::ArrayXd &dFxdy) {
JasonPries's avatar
JasonPries committed
117 118

        r = -f;
JasonPries's avatar
JasonPries committed
119
        J.setZero();
120
        for (size_t i = 0; i != TriangleQuadrature<QuadratureOrder>::size; ++i) {
JasonPries's avatar
JasonPries committed
121
            // Caclulate flux density
JasonPries's avatar
JasonPries committed
122 123
            Fx = Derivatives.dy(i).transpose() * v;
            Fy = -Derivatives.dx(i).transpose() * v;
JasonPries's avatar
JasonPries committed
124 125

            // Calculate polarization
JasonPries's avatar
JasonPries committed
126 127
            for (auto &dr : Domain.regions()) {
                dr->material().h_from_b(dr->elements(), Fx, Fy, dFxdx, dFydy, dFxdy);
JasonPries's avatar
JasonPries committed
128 129 130
            }

            // Calculate residual
JasonPries's avatar
JasonPries committed
131 132
            r += Derivatives.dy(i) * (ElementWeights(i) * Fx).matrix()
                 - Derivatives.dx(i) * (ElementWeights(i) * Fy).matrix(); // note: weak form introduces a negative sign into double-curl operator
JasonPries's avatar
JasonPries committed
133 134

            // Calculate jacobian
JasonPries's avatar
JasonPries committed
135 136 137 138 139
            J += Derivatives.dx(i) * (ElementWeights(i) * dFydy).matrix().asDiagonal() * Derivatives.dx(i).transpose()
                 + Derivatives.dy(i) * (ElementWeights(i) * dFxdx).matrix().asDiagonal() * Derivatives.dy(i).transpose()
                 + Derivatives.dx(i) * (ElementWeights(i) * dFxdy).matrix().asDiagonal() * Derivatives.dy(i).transpose()
                 + Derivatives.dy(i) * (ElementWeights(i) * dFxdy).matrix().asDiagonal() * Derivatives.dx(i).transpose();

140 141 142 143 144 145 146 147
        }
    };

    void build_matrices() {
        ElementWeights = Domain.template determinant<QuadratureOrder>();
        Basis = Domain.template basis<QuadratureOrder>();
        Derivatives = Domain.template derivative<QuadratureOrder>();

148 149
        for (size_t i = 0; i != TriangleQuadrature<QuadratureOrder>::size; ++i) {
            ElementWeights(i) *= TriangleQuadrature<QuadratureOrder>::w[i];
150 151 152
        }
    };

JasonPries's avatar
JasonPries committed
153 154
    void calculate_forcing(Eigen::VectorXd &f) {
        f.setZero();
JasonPries's avatar
JasonPries committed
155
        for (size_t i = 0; i != ForcingCondtions.size(); ++i) {
JasonPries's avatar
JasonPries committed
156 157 158 159
            f += ForcingCondtions[i]->operator()(Time);
        }
    }

JasonPries's avatar
JasonPries committed
160
    void add_current_density(std::function<double(double)> func, std::vector<std::shared_ptr<Region<2>>> regions) {
JasonPries's avatar
JasonPries committed
161
        ForcingCondtions.push_back(std::make_shared<HomogeneousForcing<2, ElementOrder, QuadratureOrder>>(func, regions, Basis, ElementWeights));
JasonPries's avatar
JasonPries committed
162 163 164 165 166 167
    }

    void remove_current_density(size_t i) {
        ForcingCondtions.erase(ForcingCondtions.begin() + i);
    }

JasonPries's avatar
JasonPries committed
168
    void add_magnetic_insulation(std::vector<std::shared_ptr<Boundary<2>>> boundaries) { // TODO: should boundaries be unique_ptrs?
JasonPries's avatar
JasonPries committed
169
        BoundaryConditions.push_back(std::make_shared<ZeroDirichlet<2, ElementOrder, QuadratureOrder>>(boundaries));
JasonPries's avatar
JasonPries committed
170 171
    }

JasonPries's avatar
JasonPries committed
172
    void add_periodic_boundary(std::vector<BoundaryMap> map, bool antiperiodic) {
JasonPries's avatar
JasonPries committed
173 174 175
        BoundaryConditions.push_back(std::make_shared<PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder>>(map, antiperiodic));
    }

176 177 178 179
    void add_periodic_boundary(std::vector<PeriodicBoundaryPair> periodic_boundary_pairs, bool antiperiodic) {
        BoundaryConditions.push_back(std::make_shared<PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder>>(periodic_boundary_pairs, Domain, antiperiodic));
    }

JasonPries's avatar
JasonPries committed
180 181 182 183 184
    template <typename... Args>
    auto add_sliding_interface(Args&&... args) {
        auto condition = std::make_shared<SlidingInterface<2, ElementOrder, QuadratureOrder>>(std::forward<Args>(args)...);
        BoundaryConditions.push_back(condition);
        return condition;
JasonPries's avatar
JasonPries committed
185 186
    }

JasonPries's avatar
JasonPries committed
187
    void apply_conditions() {
JasonPries's avatar
JasonPries committed
188
        // TODO: This is a general method that should be in a superclass
JasonPries's avatar
JasonPries committed
189 190 191
        std::vector<Eigen::Triplet<double>> triplets;
        triplets.reserve(Domain.nodes().size());

JasonPries's avatar
JasonPries committed
192 193
        for (size_t i = 0; i != Domain.nodes().size(); ++i) {
            triplets.push_back(Eigen::Triplet<double>(i, i, 1.0));
JasonPries's avatar
JasonPries committed
194 195 196 197 198 199
        }

        for (size_t i = 0; i != BoundaryConditions.size(); ++i) {
            BoundaryConditions[i]->apply(triplets);
        }

JasonPries's avatar
JasonPries committed
200
        std::set<size_t, std::less<size_t>> index;
JasonPries's avatar
JasonPries committed
201 202 203 204
        for (size_t i = 0; i != BoundaryConditions.size(); ++i) {
            BoundaryConditions[i]->reduce(index);
        }

JasonPries's avatar
JasonPries committed
205 206 207 208 209 210 211 212 213 214 215 216 217 218
        // TODO: Does this handle arbitrary number of entries in the b.c. "permutation" matrix correctly?
        // TODO: BEGIN
        // Sort triplet list by rows
        auto row_lt = [](Eigen::Triplet<double> const &x, Eigen::Triplet<double> const &y) { return x.row() < y.row(); };
        std::sort(triplets.begin(), triplets.end(), row_lt);

        size_t num_erased{0};
        auto iter{triplets.begin()};
        for (size_t const &i : index) {
            // Renumber triplets based on number of rows previously erased
            while (iter != triplets.end() && iter->row() < i) {
                *iter = Eigen::Triplet<double>(iter->row() - num_erased, iter->col(), iter->value());
                ++iter;
            }
JasonPries's avatar
JasonPries committed
219

220
            // Erase all rows associated with the current index
JasonPries's avatar
JasonPries committed
221 222
            while (iter != triplets.end() && iter->row() == i) {
                iter = triplets.erase(iter);
JasonPries's avatar
JasonPries committed
223 224
            }

JasonPries's avatar
JasonPries committed
225
            ++num_erased;
JasonPries's avatar
JasonPries committed
226 227
        }

228 229 230 231 232
        while (iter != triplets.end()) { // Finish renumbering the remainder of the entries
            *iter = Eigen::Triplet<double>(iter->row() - num_erased, iter->col(), iter->value());
            ++iter;
        }

JasonPries's avatar
JasonPries committed
233 234 235
        size_t cols{Domain.nodes().size()};
        size_t rows{cols - num_erased};
        // TODO: END
JasonPries's avatar
JasonPries committed
236

JasonPries's avatar
JasonPries committed
237 238
        Eigen::SparseMatrix<double> bc_matrix(rows, cols);
        bc_matrix.setFromTriplets(triplets.begin(), triplets.end());
JasonPries's avatar
JasonPries committed
239 240 241 242 243

        Basis.transform(bc_matrix);
        Derivatives.transform(bc_matrix);
        Unknowns = rows;
    };
244

JasonPries's avatar
JasonPries committed
245
protected:
JasonPries's avatar
JasonPries committed
246 247 248 249 250 251 252
    using PhysicsData<2, ElementOrder, QuadratureOrder>::Domain;
    using PhysicsData<2, ElementOrder, QuadratureOrder>::Derivatives;
    using PhysicsData<2, ElementOrder, QuadratureOrder>::ElementWeights;
    using PhysicsData<2, ElementOrder, QuadratureOrder>::Basis;

private:
    size_t Unknowns;
JasonPries's avatar
JasonPries committed
253
    size_t Elements;
254 255 256
};

#endif //OERSTED_PHYSICS_H