MaterialProperties.h 3.78 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 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 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
class NonlinearIsotropicMagneticMaterial : public MagneticMaterialInterface {
public:
    NonlinearIsotropicMagneticMaterial() {
        dB = 1.0;
        BSaturation = 2.0;

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

        BHSpline.push_back(std::array<double_t, 3>{a0,b0,c0});

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

        BHSpline.push_back(std::array<double_t, 3>{a1,b1,c1});

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

    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);
            double_t B = sqrt(Bx * Bx + By * By);

            double_t M, dMdB;

            if (B > BSaturation) {
                M = MSaturation;
                dMdB = 0.0;
            } else {
                size_t j = (size_t)(B / dB);
                const auto &s = BHSpline[j];

                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:
    std::vector<std::array<double_t,3>> BHSpline;
    double_t dB;
    double_t MSaturation;
    double_t BSaturation;
};

125
inline MaterialProperties Air() {
JasonPries's avatar
JasonPries committed
126 127 128 129
    return MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1.0));
}

#endif //OERSTED_MATERIALPROPERTIES_H