Physics.h 13.4 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"
13
#include "Solution.h"
JasonPries's avatar
JasonPries committed
14 15 16
#include "Node.h"
#include "PhysicsConstants.h"
#include "Triangle.h"
17

JasonPries's avatar
JasonPries committed
18 19 20 21
// 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

22
enum class FieldVariable {
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
    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
41 42 43 44 45 46 47 48 49
class PhysicsInterface {
public:
    double time() { return Time; };

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

protected:
    double Time;

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

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

JasonPries's avatar
JasonPries committed
59
    /* TODO: template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder> class PhysicsData : public PhysicsInterface {...};
JasonPries's avatar
JasonPries committed
60 61 62 63 64 65 66
    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
67
    PhysicsData(FiniteElementMesh<2, ElementOrder> &d) : Domain{d},
JasonPries's avatar
JasonPries committed
68 69 70
                                                         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
71

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

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

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

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

80 81
    Eigen::VectorXd recover_boundary(Eigen::VectorXd const &v) const { return BoundaryConditionMatrix.transpose() * v; };

JasonPries's avatar
JasonPries committed
82 83 84 85 86 87
protected:
    FiniteElementMesh<Dimension, ElementOrder> &Domain;

    DiagonalMatrixGroup<QuadraturePoints> ElementWeights;
    SparseMatrixGroup<QuadraturePoints> Basis;
    DerivativeMatrixGroup<QuadraturePoints> Derivatives;
88 89

    Eigen::SparseMatrix<double> BoundaryConditionMatrix;
90 91
};

92
template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder, FieldVariable V>
JasonPries's avatar
JasonPries committed
93
class Magnetostatic : public PhysicsData<Dimension, ElementOrder, QuadratureOrder>, public PhysicsInterface {
JasonPries's avatar
JasonPries committed
94 95 96
};

template<size_t ElementOrder, size_t QuadratureOrder>
97
class Magnetostatic<2, ElementOrder, QuadratureOrder, FieldVariable::MagneticVectorPotential> : public PhysicsData<2, ElementOrder, QuadratureOrder>, public PhysicsInterface {
98
public:
JasonPries's avatar
JasonPries committed
99
    using PhysicsData<2, ElementOrder, QuadratureOrder>::QuadraturePoints;
100

JasonPries's avatar
JasonPries committed
101 102 103 104
    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
105

106 107
    using PhysicsData<2, ElementOrder, QuadratureOrder>::recover_boundary;

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

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

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

114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139
    std::shared_ptr<Solution> initialize() {
        auto s = std::make_shared<Solution>();

        s->J = init_unknown_matrix(); // Jacobian

        s->v = init_unknown_vector(); // solution
        s->v.setZero();

        s->r = init_unknown_vector(); // residual
        s->f = init_unknown_vector(); // forcing

        s->Fx = init_element_array();
        s->Fy = init_element_array();

        s->dFxdGx = init_element_array();
        s->dFydGy = init_element_array();
        s->dFxdGy = init_element_array();

        return s;
    }

    void assemble() {
        build_matrices();
        apply_conditions();
    }

JasonPries's avatar
JasonPries committed
140 141 142
    Eigen::VectorXd init_unknown_vector() const { return Eigen::VectorXd(Unknowns); };

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

JasonPries's avatar
JasonPries committed
144 145 146 147
    Eigen::ArrayXd init_unknown_array() const { return Eigen::ArrayXd(Unknowns); };

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

148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173
    void solve(std::shared_ptr<Solution> s) {
        // TODO: Template on solver type
        // TODO: Save J,v,r,f,Fx,Fy,dFxdx,dFydy,dFxdy to reduce number of memory allocations?
        // TODO: This method is fairly general

        calculate_forcing(s->f);

        linearize(s->J, s->v, s->r, s->f, s->Fx, s->Fy, s->dFxdGx, s->dFydGy, s->dFxdGy);

        for (size_t iter = 0; iter != 20; ++iter) {
            s->Solver.compute(s->J);
            if (s->Solver.info() != Eigen::Success) {
                std::cerr << "Factorization of Jacobian failed";
            }

            s->v -= 1.0 * s->Solver.solve(s->r);
            if (s->Solver.info() != Eigen::Success) {
                std::cerr << "Solution of linear system failed";
            }

            linearize(s->J, s->v, s->r, s->f, s->Fx, s->Fy, s->dFxdGx, s->dFydGy, s->dFxdGy);

            std::cout << s->r.norm() << std::endl;
        }
    }

JasonPries's avatar
JasonPries committed
174
    void linearize(Eigen::SparseMatrix<double> &J, Eigen::VectorXd &v, Eigen::VectorXd &r, Eigen::VectorXd &f,
175
                   Eigen::ArrayXd &Fx, Eigen::ArrayXd &Fy, Eigen::ArrayXd &dFxGx, Eigen::ArrayXd &dFyGy, Eigen::ArrayXd &dFxGy) {
JasonPries's avatar
JasonPries committed
176 177

        r = -f;
JasonPries's avatar
JasonPries committed
178
        J.setZero();
179
        for (size_t i = 0; i != TriangleQuadrature<QuadratureOrder>::size; ++i) {
JasonPries's avatar
JasonPries committed
180
            // Caclulate flux density
JasonPries's avatar
JasonPries committed
181 182
            Fx = Derivatives.dy(i).transpose() * v;
            Fy = -Derivatives.dx(i).transpose() * v;
JasonPries's avatar
JasonPries committed
183

184
            // Calculate field intensity
JasonPries's avatar
JasonPries committed
185
            for (auto &dr : Domain.regions()) {
186
                dr->material().h_from_b(dr->elements(), Fx, Fy, dFxGx, dFyGy, dFxGy);
JasonPries's avatar
JasonPries committed
187 188 189
            }

            // Calculate residual
JasonPries's avatar
JasonPries committed
190 191
            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
192 193

            // Calculate jacobian
194 195 196 197
            J += Derivatives.dx(i) * (ElementWeights(i) * dFyGy).matrix().asDiagonal() * Derivatives.dx(i).transpose()
                 + Derivatives.dy(i) * (ElementWeights(i) * dFxGx).matrix().asDiagonal() * Derivatives.dy(i).transpose()
                 - Derivatives.dx(i) * (ElementWeights(i) * dFxGy).matrix().asDiagonal() * Derivatives.dy(i).transpose()
                 - Derivatives.dy(i) * (ElementWeights(i) * dFxGy).matrix().asDiagonal() * Derivatives.dx(i).transpose();
JasonPries's avatar
JasonPries committed
198

199 200 201 202 203 204 205 206
        }
    };

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

207 208
        for (size_t i = 0; i != TriangleQuadrature<QuadratureOrder>::size; ++i) {
            ElementWeights(i) *= TriangleQuadrature<QuadratureOrder>::w[i];
209 210 211
        }
    };

JasonPries's avatar
JasonPries committed
212 213
    void calculate_forcing(Eigen::VectorXd &f) {
        f.setZero();
JasonPries's avatar
JasonPries committed
214
        for (size_t i = 0; i != ForcingCondtions.size(); ++i) {
JasonPries's avatar
JasonPries committed
215 216 217 218
            f += ForcingCondtions[i]->operator()(Time);
        }
    }

219 220 221 222 223 224 225 226 227 228
    void add_current_density(std::function<double(double)> func, std::vector<std::shared_ptr<Contour const>> contours) {
        std::vector<std::shared_ptr<DiscreteRegion<2>>> regions;
        for (auto &c : contours) {
            regions.push_back(Domain.region(c));
        }

        ForcingCondtions.push_back(std::make_shared<HomogeneousForcing<2, ElementOrder, QuadratureOrder>>(func, regions, Basis, ElementWeights));
    }

    void add_current_density(std::function<double(double)> func, std::vector<std::shared_ptr<DiscreteRegion<2>>> regions) {
JasonPries's avatar
JasonPries committed
229
        ForcingCondtions.push_back(std::make_shared<HomogeneousForcing<2, ElementOrder, QuadratureOrder>>(func, regions, Basis, ElementWeights));
JasonPries's avatar
JasonPries committed
230 231 232 233 234 235
    }

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

JasonPries's avatar
JasonPries committed
236
    void add_magnetic_insulation(std::vector<std::shared_ptr<DiscreteBoundary<2>>> boundaries) { // TODO: should boundaries be unique_ptrs?
JasonPries's avatar
JasonPries committed
237
        BoundaryConditions.push_back(std::make_shared<ZeroDirichlet<2, ElementOrder, QuadratureOrder>>(boundaries));
JasonPries's avatar
JasonPries committed
238 239
    }

240 241 242 243
    void add_magnetic_insulation(std::vector<std::shared_ptr<Curve const>> boundaries) {
        BoundaryConditions.push_back(std::make_shared<ZeroDirichlet<2, ElementOrder, QuadratureOrder>>(boundaries, Domain));
    }

244 245 246 247 248 249 250 251 252 253 254
    // TODO: Using forwarding to support many constructors
    auto add_periodic_boundary(std::vector<VariableMap> map, bool antiperiodic) {
        auto pbc = std::make_shared<PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder>>(map, antiperiodic);
        BoundaryConditions.push_back(pbc);
        return pbc;
    }

    auto add_periodic_boundary(std::vector<PeriodicBoundaryPair> periodic_boundary_pairs, bool antiperiodic) {
        auto pbc = std::make_shared<PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder>>(periodic_boundary_pairs, Domain, antiperiodic);
        BoundaryConditions.push_back(pbc);
        return pbc;
JasonPries's avatar
JasonPries committed
255 256
    }

257 258 259 260
    auto add_periodic_boundary(std::vector<std::shared_ptr<DiscreteBoundary<2>>> b0, std::vector<std::shared_ptr<DiscreteBoundary<2>>> b1, std::vector<bool> orientation, bool antiperiodic) {
        auto pbc = std::make_shared<PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder>>(b0, b1, orientation, antiperiodic);
        BoundaryConditions.push_back(pbc);
        return pbc;
261 262
    }

263 264 265 266 267 268 269 270 271 272
    auto add_sliding_interface(std::shared_ptr<Curve const> interface, bool antiperiodic) {
        auto boundary = Domain.make_discontinuous(interface);
        auto condition = std::make_shared<SlidingInterface<2, ElementOrder, QuadratureOrder>>(boundary[0], boundary[1], antiperiodic);
        BoundaryConditions.push_back(condition);

        return condition;
    }

    auto add_sliding_interface(std::shared_ptr<DiscreteBoundary<2>> b0, std::shared_ptr<DiscreteBoundary<2>> b1, bool antiperiodic) {
        auto condition = std::make_shared<SlidingInterface<2, ElementOrder, QuadratureOrder>>(b0, b1, antiperiodic);
JasonPries's avatar
JasonPries committed
273 274
        BoundaryConditions.push_back(condition);
        return condition;
JasonPries's avatar
JasonPries committed
275 276
    }

JasonPries's avatar
JasonPries committed
277
    void apply_conditions() {
JasonPries's avatar
JasonPries committed
278
        // TODO: This is a general method that should be in a superclass
JasonPries's avatar
JasonPries committed
279

280 281 282
        // Apply permutation matrices
        BoundaryConditionMatrix.resize(Domain.nodes().size(), Domain.nodes().size());
        BoundaryConditionMatrix.setIdentity();
JasonPries's avatar
JasonPries committed
283 284

        for (size_t i = 0; i != BoundaryConditions.size(); ++i) {
285
            BoundaryConditions[i]->apply(BoundaryConditionMatrix);
JasonPries's avatar
JasonPries committed
286 287
        }

288 289
        // Construct reduction matrix
        std::set<size_t, std::less<size_t>> reduction_index;
JasonPries's avatar
JasonPries committed
290
        for (size_t i = 0; i != BoundaryConditions.size(); ++i) {
291
            BoundaryConditions[i]->reduce(reduction_index);
JasonPries's avatar
JasonPries committed
292 293
        }

294 295 296 297
        std::vector<Eigen::Triplet<double>> triplets;
        triplets.reserve(Domain.nodes().size());
        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
298 299
        }

300 301 302 303 304
        for (auto i = reduction_index.rbegin(); i != reduction_index.rend(); ++i) {
            triplets.erase(triplets.begin() + *i);
            for (auto j = triplets.begin() + *i; j != triplets.end(); ++j){
                *j = Eigen::Triplet<double>(j->row() - 1, j->col(), j->value()); // TODO: Could be more efficient
            }
305
        }
306 307
        Eigen::SparseMatrix<double> reduction_matrix(Domain.nodes().size() - reduction_index.size(), Domain.nodes().size());
        reduction_matrix.setFromTriplets(triplets.begin(), triplets.end());
308

309 310
        // Reduce boundary conditions
        BoundaryConditionMatrix = reduction_matrix * BoundaryConditionMatrix;
JasonPries's avatar
JasonPries committed
311

312 313
        Basis.transform(BoundaryConditionMatrix);
        Derivatives.transform(BoundaryConditionMatrix);
314
        Unknowns = Domain.nodes().size() - reduction_index.size();
JasonPries's avatar
JasonPries committed
315
    };
316

317 318 319 320
    size_t unknowns() const { return Unknowns; };

    size_t elements() const { return Elements; };

JasonPries's avatar
JasonPries committed
321
protected:
JasonPries's avatar
JasonPries committed
322 323 324 325
    using PhysicsData<2, ElementOrder, QuadratureOrder>::Domain;
    using PhysicsData<2, ElementOrder, QuadratureOrder>::Derivatives;
    using PhysicsData<2, ElementOrder, QuadratureOrder>::ElementWeights;
    using PhysicsData<2, ElementOrder, QuadratureOrder>::Basis;
326
    using PhysicsData<2, ElementOrder, QuadratureOrder>::BoundaryConditionMatrix;
JasonPries's avatar
JasonPries committed
327 328 329

private:
    size_t Unknowns;
JasonPries's avatar
JasonPries committed
330
    size_t Elements;
331 332 333
};

#endif //OERSTED_PHYSICS_H