MaterialProperties.h 4.11 KB
Newer Older
JasonPries's avatar
JasonPries committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
#ifndef OERSTED_MATERIALPROPERTIES_H
#define OERSTED_MATERIALPROPERTIES_H

#include <memory>
#include <utility>
#include <vector>

#include "Eigen"

#include "PhysicsConstants.h"

// TODO: Move into some sort of header
using Vector = Eigen::VectorXd;
using Array = Eigen::ArrayXd;

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

    virtual void h_from_b(std::vector<size_t> const &index, Array &Fx, Array &Fy, Array &dFxdx, Array &dFydy, Array &dFxdy) = 0;
};

class LinearIsotropicMagneticMaterial : public MagneticMaterialInterface {
public:
    LinearIsotropicMagneticMaterial(double relative_permeability) : Reluctivity(1.0 / (mu_0 * relative_permeability)) {};

    void h_from_b(std::vector<size_t> const &index, Array &Fx, Array &Fy, Array &dFxdx, Array &dFydy, Array &dFxdy) override {
        for(size_t const &i : index) {
            Fx(i) *= Reluctivity;
            Fy(i) *= Reluctivity;
            dFxdx(i) = Reluctivity;
            dFydy(i) = Reluctivity;
            dFxdy(i) = 0.0;
        }
    }

protected:
    double Reluctivity;
};

class MaterialProperties {
public:
    // TODO: Add std::string Name; property + a constructor which loads data saved in a file on disk
    MaterialProperties(std::shared_ptr<MagneticMaterialInterface> magnetic) : MagneticProperties{magnetic} {};

    template<typename... Args>
    void h_from_b(Args&&... args) { return MagneticProperties->h_from_b(std::forward<Args>(args)...); }

protected:
    std::shared_ptr<MagneticMaterialInterface> MagneticProperties;
};

53 54 55
class NonlinearIsotropicMagneticMaterial : public MagneticMaterialInterface {
public:
    NonlinearIsotropicMagneticMaterial() {
JasonPries's avatar
JasonPries committed
56
        DeltaB = 1.0;
57 58 59 60 61 62
        BSaturation = 2.0;

        double_t a0 = 0.0;
        double_t b0 = 0.999 / mu_0;
        double_t c0 = 0.0;

JasonPries's avatar
JasonPries committed
63
        MBSpline.push_back(std::array<double_t, 3>{a0,b0,c0});
64 65 66 67 68

        double_t a1 = -b0 / 2.0;
        double_t b1 = -4.0 * a1;
        double_t c1 = b0 - a1 - b1;

JasonPries's avatar
JasonPries committed
69
        MBSpline.push_back(std::array<double_t, 3>{a1,b1,c1});
70 71 72 73

        MSaturation = (a1 * BSaturation + b1) * BSaturation + c1;
    };

JasonPries's avatar
JasonPries committed
74 75 76 77 78
    NonlinearIsotropicMagneticMaterial(std::vector<std::array<double_t,3>> mbspline, double_t db, double_t bsat) : MBSpline(mbspline), DeltaB(db), BSaturation(bsat) {
        const auto &s = MBSpline.back();
        MSaturation = (s[0] * BSaturation + s[1]) * BSaturation + s[2];
    }

79 80 81 82
    void h_from_b(std::vector<size_t> const &index, Array &Fx, Array &Fy, Array &dFxdx, Array &dFydy, Array &dFxdy) override {
        for(size_t const &i : index) {
            const double_t &Bx = Fx(i);
            const double_t &By = Fy(i);
JasonPries's avatar
JasonPries committed
83

84 85 86 87 88 89 90 91
            double_t B = sqrt(Bx * Bx + By * By);

            double_t M, dMdB;

            if (B > BSaturation) {
                M = MSaturation;
                dMdB = 0.0;
            } else {
JasonPries's avatar
JasonPries committed
92 93
                size_t j = (size_t)(B / DeltaB);
                const auto &s = MBSpline[j];
94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124

                M = (s[0] * B + s[1]) * B + s[2];
                dMdB = (2.0 * s[0] * B + s[1]);
            }

            if (B > 0) {
                double_t Mx = Bx * (M / B);
                double_t My = By * (M / B);

                double_t dMxdBx = M / B + (Bx / B) * (Bx / B) * (dMdB - M / B);
                double_t dMydBy = M / B + (By / B) * (By / B) * (dMdB - M / B);
                double_t dMxdBy = (Bx / B) * (By / B) * (dMdB - M / B);

                dFxdx(i) = 1.0 / mu_0 - dMxdBx;
                dFydy(i) = 1.0 / mu_0 - dMydBy;
                dFxdy(i) = -dMxdBy;

                Fx(i) = Bx / mu_0 - Mx;
                Fy(i) = By / mu_0 - My;
            } else {
                dFxdx(i) = 1.0 / mu_0 - dMdB;
                dFydy(i) = 1.0 / mu_0 - dMdB;
                dFxdy(i) = 0.0;

                Fx(i) = 0.0;
                Fy(i) = 0.0;
            }
        }
    }

protected:
JasonPries's avatar
JasonPries committed
125 126
    std::vector<std::array<double_t,3>> MBSpline;
    double_t DeltaB;
127 128 129 130
    double_t MSaturation;
    double_t BSaturation;
};

131
inline MaterialProperties Air() {
JasonPries's avatar
JasonPries committed
132 133 134
    return MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1.0));
}

JasonPries's avatar
JasonPries committed
135 136
MaterialProperties JFE_35JN210();

JasonPries's avatar
JasonPries committed
137
#endif //OERSTED_MATERIALPROPERTIES_H