Physics.h 13.6 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
    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);

JasonPries's avatar
JasonPries committed
155

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

JasonPries's avatar
JasonPries committed
158 159 160 161
        size_t iter = 0;
        double_t norm_dv = 1.0;

        while (iter < 30 && norm_dv > sqrt(DBL_EPSILON) * std::sqrt(s->v.size())) {
162 163 164 165 166
            s->Solver.compute(s->J);
            if (s->Solver.info() != Eigen::Success) {
                std::cerr << "Factorization of Jacobian failed";
            }

JasonPries's avatar
JasonPries committed
167
            auto dv = s->Solver.solve(s->r);
168 169 170
            if (s->Solver.info() != Eigen::Success) {
                std::cerr << "Solution of linear system failed";
            }
JasonPries's avatar
JasonPries committed
171 172 173
            s->v -= dv;

            norm_dv = dv.norm() / s->v.norm();
174 175 176

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

JasonPries's avatar
JasonPries committed
177
            iter++;
178 179 180
        }
    }

JasonPries's avatar
JasonPries committed
181
    void linearize(Eigen::SparseMatrix<double> &J, Eigen::VectorXd &v, Eigen::VectorXd &r, Eigen::VectorXd &f,
182
                   Eigen::ArrayXd &Fx, Eigen::ArrayXd &Fy, Eigen::ArrayXd &dFxGx, Eigen::ArrayXd &dFyGy, Eigen::ArrayXd &dFxGy) {
JasonPries's avatar
JasonPries committed
183 184

        r = -f;
JasonPries's avatar
JasonPries committed
185
        J.setZero();
186
        for (size_t i = 0; i != TriangleQuadrature<QuadratureOrder>::size; ++i) {
JasonPries's avatar
JasonPries committed
187
            // Caclulate flux density
JasonPries's avatar
JasonPries committed
188 189
            Fx = Derivatives.dy(i).transpose() * v;
            Fy = -Derivatives.dx(i).transpose() * v;
JasonPries's avatar
JasonPries committed
190

191
            // Calculate field intensity
JasonPries's avatar
JasonPries committed
192
            for (auto &dr : Domain.regions()) {
193
                dr->material().h_from_b(dr->elements(), Fx, Fy, dFxGx, dFyGy, dFxGy);
JasonPries's avatar
JasonPries committed
194 195 196
            }

            // Calculate residual
JasonPries's avatar
JasonPries committed
197 198
            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
199 200

            // Calculate jacobian
201 202 203 204
            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
205

206 207 208 209 210 211 212 213
        }
    };

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

214 215
        for (size_t i = 0; i != TriangleQuadrature<QuadratureOrder>::size; ++i) {
            ElementWeights(i) *= TriangleQuadrature<QuadratureOrder>::w[i];
216 217 218
        }
    };

JasonPries's avatar
JasonPries committed
219 220
    void calculate_forcing(Eigen::VectorXd &f) {
        f.setZero();
JasonPries's avatar
JasonPries committed
221
        for (size_t i = 0; i != ForcingCondtions.size(); ++i) {
JasonPries's avatar
JasonPries committed
222 223 224 225
            f += ForcingCondtions[i]->operator()(Time);
        }
    }

226 227 228 229 230 231 232 233 234 235
    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
236
        ForcingCondtions.push_back(std::make_shared<HomogeneousForcing<2, ElementOrder, QuadratureOrder>>(func, regions, Basis, ElementWeights));
JasonPries's avatar
JasonPries committed
237 238 239 240 241 242
    }

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

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

247 248 249 250
    void add_magnetic_insulation(std::vector<std::shared_ptr<Curve const>> boundaries) {
        BoundaryConditions.push_back(std::make_shared<ZeroDirichlet<2, ElementOrder, QuadratureOrder>>(boundaries, Domain));
    }

251 252 253 254 255 256 257 258 259 260 261
    // 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
262 263
    }

264 265 266 267
    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;
268 269
    }

270 271 272 273 274 275 276 277 278 279
    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
280 281
        BoundaryConditions.push_back(condition);
        return condition;
JasonPries's avatar
JasonPries committed
282 283
    }

JasonPries's avatar
JasonPries committed
284
    void apply_conditions() {
JasonPries's avatar
JasonPries committed
285
        // TODO: This is a general method that should be in a superclass
JasonPries's avatar
JasonPries committed
286

287 288 289
        // Apply permutation matrices
        BoundaryConditionMatrix.resize(Domain.nodes().size(), Domain.nodes().size());
        BoundaryConditionMatrix.setIdentity();
JasonPries's avatar
JasonPries committed
290 291

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

295 296
        // Construct reduction matrix
        std::set<size_t, std::less<size_t>> reduction_index;
JasonPries's avatar
JasonPries committed
297
        for (size_t i = 0; i != BoundaryConditions.size(); ++i) {
298
            BoundaryConditions[i]->reduce(reduction_index);
JasonPries's avatar
JasonPries committed
299 300
        }

301 302 303 304
        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
305 306
        }

307 308 309 310 311
        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
            }
312
        }
313 314
        Eigen::SparseMatrix<double> reduction_matrix(Domain.nodes().size() - reduction_index.size(), Domain.nodes().size());
        reduction_matrix.setFromTriplets(triplets.begin(), triplets.end());
315

316 317
        // Reduce boundary conditions
        BoundaryConditionMatrix = reduction_matrix * BoundaryConditionMatrix;
JasonPries's avatar
JasonPries committed
318

319 320
        Basis.transform(BoundaryConditionMatrix);
        Derivatives.transform(BoundaryConditionMatrix);
321
        Unknowns = Domain.nodes().size() - reduction_index.size();
JasonPries's avatar
JasonPries committed
322
    };
323

324 325 326 327
    size_t unknowns() const { return Unknowns; };

    size_t elements() const { return Elements; };

JasonPries's avatar
JasonPries committed
328
protected:
JasonPries's avatar
JasonPries committed
329 330 331 332
    using PhysicsData<2, ElementOrder, QuadratureOrder>::Domain;
    using PhysicsData<2, ElementOrder, QuadratureOrder>::Derivatives;
    using PhysicsData<2, ElementOrder, QuadratureOrder>::ElementWeights;
    using PhysicsData<2, ElementOrder, QuadratureOrder>::Basis;
333
    using PhysicsData<2, ElementOrder, QuadratureOrder>::BoundaryConditionMatrix;
JasonPries's avatar
JasonPries committed
334 335 336

private:
    size_t Unknowns;
JasonPries's avatar
JasonPries committed
337
    size_t Elements;
338 339 340
};

#endif //OERSTED_PHYSICS_H