test_FiniteElementMesh.cpp 38.5 KB
Newer Older
1 2
#include "test_Physics.hpp"

3 4 5
class MasterTriangle : public ::testing::Test {
public:
    virtual void SetUp() {
JasonPries's avatar
JasonPries committed
6 7 8
        node.push_back(XY{0.0, 0.0});
        node.push_back(XY{1.0, 0.0});
        node.push_back(XY{0.0, 1.0});
9

10
        tri.push_back(Triangle<1>(0, {1, 2, 0}));
11

JasonPries's avatar
JasonPries committed
12
        reg.push_back(std::make_shared<Region<2>>(std::vector<size_t>{0}));
JasonPries's avatar
JasonPries committed
13

JasonPries's avatar
JasonPries committed
14 15 16
        bnd.push_back(std::make_shared<Boundary<2>>(std::vector<size_t>{0, 1}));
        bnd.push_back(std::make_shared<Boundary<2>>(std::vector<size_t>{0, 2}));
        bnd.push_back(std::make_shared<Boundary<2>>(std::vector<size_t>{1, 2}));
JasonPries's avatar
JasonPries committed
17

JasonPries's avatar
JasonPries committed
18
        femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
19 20 21 22
    }

    std::vector<XY> node;
    std::vector<Triangle<1>> tri;
JasonPries's avatar
JasonPries committed
23 24 25

    std::vector<std::shared_ptr<Region<2>>> reg;
    std::vector<std::shared_ptr<Boundary<2>>> bnd;
26

JasonPries's avatar
JasonPries committed
27
    FiniteElementMesh<2, 1> femesh;
28 29 30 31 32 33
};

TEST_F(MasterTriangle, mass_matrix) {
    auto mat = femesh.basis<2>();
    auto det = femesh.determinant<2>();

34 35 36
    Eigen::SparseMatrix<double> M = (mat[0] * det(0).matrix().asDiagonal() * mat[0].transpose()) * TriangleQuadrature<2>::w[0]
                                    + (mat[1] * det(1).matrix().asDiagonal() * mat[1].transpose()) * TriangleQuadrature<2>::w[1]
                                    + (mat[2] * det(2).matrix().asDiagonal() * mat[2].transpose()) * TriangleQuadrature<2>::w[2];
37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52

    EXPECT_NEAR(M.coeffRef(0, 0), 1.0 / 12.0, FLT_EPSILON);
    EXPECT_NEAR(M.coeffRef(1, 0), 1.0 / 24.0, FLT_EPSILON);
    EXPECT_NEAR(M.coeffRef(2, 0), 1.0 / 24.0, FLT_EPSILON);
    EXPECT_NEAR(M.coeffRef(0, 1), 1.0 / 24.0, FLT_EPSILON);
    EXPECT_NEAR(M.coeffRef(1, 1), 1.0 / 12.0, FLT_EPSILON);
    EXPECT_NEAR(M.coeffRef(2, 1), 1.0 / 24.0, FLT_EPSILON);
    EXPECT_NEAR(M.coeffRef(0, 2), 1.0 / 24.0, FLT_EPSILON);
    EXPECT_NEAR(M.coeffRef(1, 2), 1.0 / 24.0, FLT_EPSILON);
    EXPECT_NEAR(M.coeffRef(2, 2), 1.0 / 12.0, FLT_EPSILON);
}

TEST_F(MasterTriangle, stiffness_matrix) {
    auto df = femesh.derivative<1>();
    auto det = femesh.determinant<1>();

JasonPries's avatar
JasonPries committed
53
    Eigen::SparseMatrix<double> K = (df.dx(0) * det(0).matrix().asDiagonal() * df.dx(0).transpose()
54
                                     + df.dy(0) * det(0).matrix().asDiagonal() * df.dy(0).transpose()) * TriangleQuadrature<1>::w[0];
55 56 57 58 59 60 61 62 63 64 65 66 67

    EXPECT_NEAR(K.coeffRef(0, 0), 1.0, FLT_EPSILON);
    EXPECT_NEAR(K.coeffRef(1, 0), -1.0 / 2.0, FLT_EPSILON);
    EXPECT_NEAR(K.coeffRef(2, 0), -1.0 / 2.0, FLT_EPSILON);
    EXPECT_NEAR(K.coeffRef(0, 1), -1.0 / 2.0, FLT_EPSILON);
    EXPECT_NEAR(K.coeffRef(1, 1), 1.0 / 2.0, FLT_EPSILON);
    EXPECT_NEAR(K.coeffRef(2, 1), 0.0, FLT_EPSILON);
    EXPECT_NEAR(K.coeffRef(0, 2), -1.0 / 2.0, FLT_EPSILON);
    EXPECT_NEAR(K.coeffRef(1, 2), 0.0, FLT_EPSILON);
    EXPECT_NEAR(K.coeffRef(2, 2), 1.0 / 2.0, FLT_EPSILON);
}

TEST_F(MasterTriangle, magnetostatic) {
JasonPries's avatar
JasonPries committed
68
    //TODO: This isn't really testing anything...
JasonPries's avatar
JasonPries committed
69 70
    static constexpr size_t Q = 1;

71
    Magnetostatic<2, 1, Q, FieldVariable::A> ms(femesh);
JasonPries's avatar
JasonPries committed
72
    ms.time(1.0);
JasonPries's avatar
JasonPries committed
73

JasonPries's avatar
JasonPries committed
74
    auto ff = [](double t) { return 1.0 * t; };
JasonPries's avatar
JasonPries committed
75

JasonPries's avatar
JasonPries committed
76
    ms.add_current_density(ff, reg);
77 78 79

    ms.build_matrices();

JasonPries's avatar
JasonPries committed
80 81 82 83
    auto J = ms.init_unknown_matrix();
    auto v = ms.init_unknown_vector();
    auto r = ms.init_unknown_vector();
    auto f = ms.init_unknown_vector();
JasonPries's avatar
JasonPries committed
84 85 86 87 88
    auto Fx = ms.init_element_array();
    auto Fy = ms.init_element_array();
    auto dFxdx = ms.init_element_array();
    auto dFydy = ms.init_element_array();
    auto dFxdy = ms.init_element_array();
89

JasonPries's avatar
JasonPries committed
90
    ms.calculate_forcing(f);
JasonPries's avatar
JasonPries committed
91
    ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
92

JasonPries's avatar
JasonPries committed
93
    //auto ff = [](double t) { return 1.0 * t; };
94

JasonPries's avatar
JasonPries committed
95
    //HomogeneousForcing<2,1,Q> hf{ff, std::vector<size_t>{0}, ms.basis(), ms.weights(), ms.domain()};
96

JasonPries's avatar
JasonPries committed
97
    std::cout << f << std::endl;
98 99 100
}

class MasterTriangleRotationReflection : public ::testing::Test {
101 102 103 104 105 106 107 108
public:
    virtual void SetUp() {
        node.push_back(XY{0.0, 0.0});
        node.push_back(XY{1.0, 0.0});
        node.push_back(XY{0.0, 1.0});
        node.push_back(XY{-1.0, 0.0});
        node.push_back(XY{0.0, -1.0});

109 110 111 112 113 114 115
        tri.push_back(Triangle<1>(0, {1, 2, 0})); // Master
        tri.push_back(Triangle<1>(1, {2, 3, 0})); // 90 degree
        tri.push_back(Triangle<1>(2, {3, 4, 0})); // 180 degree
        tri.push_back(Triangle<1>(3, {4, 1, 0})); // 270 degree
        tri.push_back(Triangle<1>(4, {3, 2, 0})); // y-axis reflection
        tri.push_back(Triangle<1>(5, {1, 4, 0})); // x-axis reflection
        tri.push_back(Triangle<1>(6, {4, 3, 0})); // x=y reflection
JasonPries's avatar
JasonPries committed
116

JasonPries's avatar
JasonPries committed
117
        reg.push_back(std::make_shared<Region<2>>(std::vector<size_t>{0, 1, 2, 3})); // Triangles 4,5,6 are duplicates
JasonPries's avatar
JasonPries committed
118

JasonPries's avatar
JasonPries committed
119
        femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
120 121 122 123
    }

    std::vector<XY> node;
    std::vector<Triangle<1>> tri;
JasonPries's avatar
JasonPries committed
124

JasonPries's avatar
JasonPries committed
125 126
    std::vector<std::shared_ptr<Region<2>>> reg;
    std::vector<std::shared_ptr<Boundary<2>>> bnd;
JasonPries's avatar
JasonPries committed
127 128

    FiniteElementMesh<2, 1> femesh;
129 130
};

131
TEST_F(MasterTriangleRotationReflection, jacobian_1) {
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 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 174 175 176 177 178 179 180 181
    auto J = tri[0].jacobian<1>(node);

    EXPECT_NEAR(J(0, 0), 1.0, FLT_EPSILON);
    EXPECT_NEAR(J(0, 1), 0.0, FLT_EPSILON);
    EXPECT_NEAR(J(1, 0), 0.0, FLT_EPSILON);
    EXPECT_NEAR(J(1, 1), 1.0, FLT_EPSILON);

    J = tri[1].jacobian<1>(node);

    EXPECT_NEAR(J(0, 0), 0.0, FLT_EPSILON);
    EXPECT_NEAR(J(0, 1), -1.0, FLT_EPSILON);
    EXPECT_NEAR(J(1, 0), 1.0, FLT_EPSILON);
    EXPECT_NEAR(J(1, 1), 0.0, FLT_EPSILON);

    J = tri[2].jacobian<1>(node);

    EXPECT_NEAR(J(0, 0), -1.0, FLT_EPSILON);
    EXPECT_NEAR(J(0, 1), 0.0, FLT_EPSILON);
    EXPECT_NEAR(J(1, 0), 0.0, FLT_EPSILON);
    EXPECT_NEAR(J(1, 1), -1.0, FLT_EPSILON);

    J = tri[3].jacobian<1>(node);

    EXPECT_NEAR(J(0, 0), 0.0, FLT_EPSILON);
    EXPECT_NEAR(J(0, 1), 1.0, FLT_EPSILON);
    EXPECT_NEAR(J(1, 0), -1.0, FLT_EPSILON);
    EXPECT_NEAR(J(1, 1), 0.0, FLT_EPSILON);

    J = tri[4].jacobian<1>(node);

    EXPECT_NEAR(J(0, 0), -1.0, FLT_EPSILON);
    EXPECT_NEAR(J(0, 1), 0.0, FLT_EPSILON);
    EXPECT_NEAR(J(1, 0), 0.0, FLT_EPSILON);
    EXPECT_NEAR(J(1, 1), 1.0, FLT_EPSILON);

    J = tri[5].jacobian<1>(node);

    EXPECT_NEAR(J(0, 0), 1.0, FLT_EPSILON);
    EXPECT_NEAR(J(0, 1), 0.0, FLT_EPSILON);
    EXPECT_NEAR(J(1, 0), 0.0, FLT_EPSILON);
    EXPECT_NEAR(J(1, 1), -1.0, FLT_EPSILON);

    J = tri[6].jacobian<1>(node);

    EXPECT_NEAR(J(0, 0), 0.0, FLT_EPSILON);
    EXPECT_NEAR(J(0, 1), -1.0, FLT_EPSILON);
    EXPECT_NEAR(J(1, 0), -1.0, FLT_EPSILON);
    EXPECT_NEAR(J(1, 1), 0.0, FLT_EPSILON);
}

182
TEST_F(MasterTriangleRotationReflection, basis_1) {
JasonPries's avatar
JasonPries committed
183
    auto mat = femesh.basis<1>();
184 185

    for (size_t i = 0; i != tri.size(); ++i) {
JasonPries's avatar
JasonPries committed
186 187 188
        EXPECT_NEAR(mat(0, tri[i].node(0), i), 1.0 / 3.0, FLT_EPSILON);
        EXPECT_NEAR(mat(0, tri[i].node(1), i), 1.0 / 3.0, FLT_EPSILON);
        EXPECT_NEAR(mat(0, tri[i].node(2), i), 1.0 / 3.0, FLT_EPSILON);
189 190 191
    }
}

192
TEST_F(MasterTriangleRotationReflection, derivative_1) {
JasonPries's avatar
JasonPries committed
193
    auto df = femesh.derivative<1>();
194

JasonPries's avatar
JasonPries committed
195 196 197 198 199 200
    EXPECT_NEAR(df.dx(0, tri[0].node(0), 0), 1.0, FLT_EPSILON);
    EXPECT_NEAR(df.dy(0, tri[0].node(0), 0), 0.0, FLT_EPSILON);
    EXPECT_NEAR(df.dx(0, tri[0].node(1), 0), 0.0, FLT_EPSILON);
    EXPECT_NEAR(df.dy(0, tri[0].node(1), 0), 1.0, FLT_EPSILON);
    EXPECT_NEAR(df.dx(0, tri[0].node(2), 0), -1.0, FLT_EPSILON);
    EXPECT_NEAR(df.dy(0, tri[0].node(2), 0), -1.0, FLT_EPSILON);
201 202
}

203
TEST_F(MasterTriangleRotationReflection, basis_2) {
JasonPries's avatar
JasonPries committed
204
    auto mat = femesh.basis<2>();
205 206

    for (size_t i = 0; i != tri.size(); ++i) {
JasonPries's avatar
JasonPries committed
207 208 209 210 211 212 213
        EXPECT_NEAR(mat(0, tri[i].node(0), i), 2.0 / 3.0, FLT_EPSILON);
        EXPECT_NEAR(mat(0, tri[i].node(1), i), 1.0 / 6.0, FLT_EPSILON);
        EXPECT_NEAR(mat(0, tri[i].node(1), i), 1.0 / 6.0, FLT_EPSILON);

        EXPECT_NEAR(mat(1, tri[i].node(0), i), 1.0 / 6.0, FLT_EPSILON);
        EXPECT_NEAR(mat(1, tri[i].node(1), i), 2.0 / 3.0, FLT_EPSILON);
        EXPECT_NEAR(mat(1, tri[i].node(2), i), 1.0 / 6.0, FLT_EPSILON);
214

JasonPries's avatar
JasonPries committed
215 216 217 218 219
        EXPECT_NEAR(mat(2, tri[i].node(0), i), 1.0 / 6.0, FLT_EPSILON);
        EXPECT_NEAR(mat(2, tri[i].node(1), i), 1.0 / 6.0, FLT_EPSILON);
        EXPECT_NEAR(mat(2, tri[i].node(2), i), 2.0 / 3.0, FLT_EPSILON);
    }
}
220

221
TEST_F(MasterTriangleRotationReflection, derivative_2) {
JasonPries's avatar
JasonPries committed
222
    auto df = femesh.derivative<2>();
223

JasonPries's avatar
JasonPries committed
224 225 226 227 228 229 230
    for (size_t q = 0; q != 3; ++q) {
        EXPECT_NEAR(df.dx(q, tri[0].node(0), 0), 1.0, FLT_EPSILON);
        EXPECT_NEAR(df.dy(q, tri[0].node(0), 0), 0.0, FLT_EPSILON);
        EXPECT_NEAR(df.dx(q, tri[0].node(1), 0), 0.0, FLT_EPSILON);
        EXPECT_NEAR(df.dy(q, tri[0].node(1), 0), 1.0, FLT_EPSILON);
        EXPECT_NEAR(df.dx(q, tri[0].node(2), 0), -1.0, FLT_EPSILON);
        EXPECT_NEAR(df.dy(q, tri[0].node(2), 0), -1.0, FLT_EPSILON);
231 232 233
    }
}

234
class SimpleSquare : public ::testing::Test {
JasonPries's avatar
JasonPries committed
235 236 237 238 239 240 241
public:
    virtual void SetUp() {
        node.push_back(XY{0.0, 0.0});
        node.push_back(XY{1.0, 0.0});
        node.push_back(XY{0.0, 1.0});
        node.push_back(XY{1.0, 1.0});

242 243
        tri.push_back(Triangle<1>(0, {0, 1, 2}));
        tri.push_back(Triangle<1>(1, {1, 3, 2}));
244

JasonPries's avatar
JasonPries committed
245 246
        reg.push_back(std::make_shared<Region<2>>(std::vector<size_t>{0}));
        reg.push_back(std::make_shared<Region<2>>(std::vector<size_t>{1}));
JasonPries's avatar
JasonPries committed
247

JasonPries's avatar
JasonPries committed
248 249 250 251 252
        bnd.push_back(std::make_shared<Boundary<2>>(std::vector<size_t>{0, 1}));
        bnd.push_back(std::make_shared<Boundary<2>>(std::vector<size_t>{0, 2}));
        bnd.push_back(std::make_shared<Boundary<2>>(std::vector<size_t>{1, 2}));
        bnd.push_back(std::make_shared<Boundary<2>>(std::vector<size_t>{1, 3}));
        bnd.push_back(std::make_shared<Boundary<2>>(std::vector<size_t>{2, 3}));
JasonPries's avatar
JasonPries committed
253 254

        femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
JasonPries's avatar
JasonPries committed
255 256 257 258 259

        vx_m << 1.0, 0.0, 1.0, 0.0;
        vx_p << 0.0, 1.0, 0.0, 1.0;
        vy_m << 1.0, 1.0, 0.0, 0.0;
        vy_p << 0.0, 0.0, 1.0, 1.0;
260 261
    }

JasonPries's avatar
JasonPries committed
262 263
    std::vector<XY> node;
    std::vector<Triangle<1>> tri;
JasonPries's avatar
JasonPries committed
264 265 266

    std::vector<std::shared_ptr<Region<2>>> reg;
    std::vector<std::shared_ptr<Boundary<2>>> bnd;
JasonPries's avatar
JasonPries committed
267 268 269 270 271 272 273 274 275

    Eigen::Vector4d vx_m;
    Eigen::Vector4d vx_p;
    Eigen::Vector4d vy_m;
    Eigen::Vector4d vy_p;

    FiniteElementMesh<2, 1> femesh;
};

276
TEST_F(SimpleSquare, derivative_1_temp) {
JasonPries's avatar
JasonPries committed
277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313
    auto df = femesh.derivative<1>();

    Eigen::Vector2d tmp;

    tmp = df.dx(0).transpose() * vx_p;
    EXPECT_NEAR(tmp(0), 1.0, FLT_EPSILON);
    EXPECT_NEAR(tmp(0), 1.0, FLT_EPSILON);

    tmp = df.dx(0).transpose() * vx_m;
    EXPECT_NEAR(tmp(0), -1.0, FLT_EPSILON);
    EXPECT_NEAR(tmp(0), -1.0, FLT_EPSILON);

    tmp = df.dx(0).transpose() * vy_m;
    EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
    EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);

    tmp = df.dx(0).transpose() * vy_m;
    EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
    EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);

    tmp = df.dy(0).transpose() * vy_p;
    EXPECT_NEAR(tmp(0), 1.0, FLT_EPSILON);
    EXPECT_NEAR(tmp(0), 1.0, FLT_EPSILON);

    tmp = df.dy(0).transpose() * vy_m;
    EXPECT_NEAR(tmp(0), -1.0, FLT_EPSILON);
    EXPECT_NEAR(tmp(0), -1.0, FLT_EPSILON);

    tmp = df.dy(0).transpose() * vx_m;
    EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
    EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);

    tmp = df.dy(0).transpose() * vx_m;
    EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
    EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
}

314
TEST_F(SimpleSquare, derivative_2_temp) {
JasonPries's avatar
JasonPries committed
315 316 317
    auto df = femesh.derivative<2>();

    Eigen::Vector2d tmp;
318 319

    for (size_t i = 0; i != 3; ++i) {
JasonPries's avatar
JasonPries committed
320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350
        tmp = df.dx(i).transpose() * vx_p;
        EXPECT_NEAR(tmp(0), 1.0, FLT_EPSILON);
        EXPECT_NEAR(tmp(0), 1.0, FLT_EPSILON);

        tmp = df.dx(i).transpose() * vx_m;
        EXPECT_NEAR(tmp(0), -1.0, FLT_EPSILON);
        EXPECT_NEAR(tmp(0), -1.0, FLT_EPSILON);

        tmp = df.dx(i).transpose() * vy_m;
        EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
        EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);

        tmp = df.dx(i).transpose() * vy_m;
        EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
        EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);

        tmp = df.dy(i).transpose() * vy_p;
        EXPECT_NEAR(tmp(0), 1.0, FLT_EPSILON);
        EXPECT_NEAR(tmp(0), 1.0, FLT_EPSILON);

        tmp = df.dy(i).transpose() * vy_m;
        EXPECT_NEAR(tmp(0), -1.0, FLT_EPSILON);
        EXPECT_NEAR(tmp(0), -1.0, FLT_EPSILON);

        tmp = df.dy(i).transpose() * vx_m;
        EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
        EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);

        tmp = df.dy(i).transpose() * vx_m;
        EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
        EXPECT_NEAR(tmp(0), 0.0, FLT_EPSILON);
351
    }
JasonPries's avatar
JasonPries committed
352 353 354 355 356 357 358 359 360 361
}

TEST_F(SimpleSquare, magnetostatic_forcing) {
    // TODO: Write Physics::add_condition() (templated?)
    // TODO: BoundaryCondition, ForcingCondition, MotionCondition?
    // TODO: ?then template parameters of conditions can be hidden?
    // TODO: Hierarchy of conditions: Those altering RHS only (forcing, some bcs), those altering matrices in a fixed way (some bcs), those altering the mesh (motion, some bcs)
    //      TODO: Applied in reverse order, Motion->BoundaryCondition->Forcing,
    //      TODO: Or, TimeVarying->TimeInvariant->Forcing?

362
    Magnetostatic<2, 1, 2, FieldVariable::A> ms{femesh};
JasonPries's avatar
JasonPries committed
363
    ms.time(1.0);
JasonPries's avatar
JasonPries committed
364 365 366
    ms.build_matrices();

    auto ff = [](double t) { return 1.0 * t; };
JasonPries's avatar
JasonPries committed
367
    auto f = ms.init_unknown_vector();
JasonPries's avatar
JasonPries committed
368

JasonPries's avatar
JasonPries committed
369
    ms.add_current_density(ff, {reg[0]});
JasonPries's avatar
JasonPries committed
370
    ms.calculate_forcing(f);
JasonPries's avatar
JasonPries committed
371 372 373 374
    EXPECT_NEAR(f(0), 1.0 / 6.0, FLT_EPSILON);
    EXPECT_NEAR(f(1), 1.0 / 6.0, FLT_EPSILON);
    EXPECT_NEAR(f(2), 1.0 / 6.0, FLT_EPSILON);
    EXPECT_NEAR(f(3), 0.0 / 6.0, FLT_EPSILON);
JasonPries's avatar
JasonPries committed
375 376
    ms.remove_current_density(0);

JasonPries's avatar
JasonPries committed
377
    ms.add_current_density(ff, {reg[1]});
JasonPries's avatar
JasonPries committed
378
    ms.calculate_forcing(f);
JasonPries's avatar
JasonPries committed
379 380 381 382
    EXPECT_NEAR(f(0), 0.0 / 6.0, FLT_EPSILON);
    EXPECT_NEAR(f(1), 1.0 / 6.0, FLT_EPSILON);
    EXPECT_NEAR(f(2), 1.0 / 6.0, FLT_EPSILON);
    EXPECT_NEAR(f(3), 1.0 / 6.0, FLT_EPSILON);
JasonPries's avatar
JasonPries committed
383 384
    ms.remove_current_density(0);

JasonPries's avatar
JasonPries committed
385
    ms.add_current_density(ff, {reg[0], reg[1]});
JasonPries's avatar
JasonPries committed
386
    ms.calculate_forcing(f);
JasonPries's avatar
JasonPries committed
387 388 389 390
    EXPECT_NEAR(f(0), 1.0 / 6.0, FLT_EPSILON);
    EXPECT_NEAR(f(1), 2.0 / 6.0, FLT_EPSILON);
    EXPECT_NEAR(f(2), 2.0 / 6.0, FLT_EPSILON);
    EXPECT_NEAR(f(3), 1.0 / 6.0, FLT_EPSILON);
JasonPries's avatar
JasonPries committed
391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449
    ms.remove_current_density(0);

    /*
    auto J = ms.init_matrix();
    auto v = ms.init_vector();
    auto r = ms.init_vector();
    v.setZero();

    ms.linearize(J,v,r,f);

    //Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
    //solver.compute(J);
    //if(solver.info()!=Eigen::Success) {
    //    std::cout << "Failed" << std::endl;
    //} else {
    //    std::cout << "Success" << std::endl;
    //}

    Eigen::SparseMatrix<double> J3 = J.block(1,1,3,3);
    Eigen::SparseMatrix<double> J2 = J.block(1,0,3,1);
    Eigen::SparseMatrix<double> J1 = J.block(0,1,1,3);
    Eigen::SparseMatrix<double> J0 = J.block(0,0,1,1);

    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
    Eigen::SparseMatrix<double> Jb = J3 - (J2 * (J0 \ J1));

    std::cout << J3 << std::endl;
    std::cout << J2 << std::endl;
    std::cout << J1 << std::endl;
    std::cout << J0 << std::endl;
    std::cout << Jb << std::endl;

    Eigen::ConjugateGradient<Eigen::SparseMatrix<double>, Eigen::Lower|Eigen::Upper> cg;
    cg.compute(J);
    if(cg.info()!=Eigen::Success) {
        std::cout << "Failed" << std::endl;
    } else {
        std::cout << "Success" << std::endl;
    }
    v -= cg.solve(r);

    std::cout << J << std::endl;
    std::cout << v << std::endl;
    std::cout << r << std::endl;
    std::cout << f << std::endl;
    std::cout << J * v - r << std::endl;


    ms.linearize(J,v,r,f);
    std::cout << J << std::endl;
    std::cout << v << std::endl;
    std::cout << r << std::endl;
    std::cout << f << std::endl;

    std::cout << ms.derivatives().dx(0).transpose() * v << std::endl;
    std::cout << ms.derivatives().dy(0).transpose() * v << std::endl;
    */
}

JasonPries's avatar
JasonPries committed
450
class TwoRegionHexagon : public ::testing::Test {
JasonPries's avatar
JasonPries committed
451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466
public:
    virtual void SetUp() {
        node.push_back(XY{0.0, 0.0});

        double r{1.0};
        for (size_t i = 0; i != 6; ++i) {
            double a = (i * M_PI) / 3.0;
            node.push_back(XY{r * std::cos(a), r * std::sin(a)});
        }

        r = 2.0;
        for (size_t i = 0; i != 6; ++i) {
            double a = (i * M_PI) / 3.0 + M_PI / 6.0;
            node.push_back(XY{r * std::cos(a), r * std::sin(a)});
        }

467 468 469 470 471 472
        tri.push_back(Triangle<1>(0, {0, 1, 2}));
        tri.push_back(Triangle<1>(1, {0, 2, 3}));
        tri.push_back(Triangle<1>(2, {0, 3, 4}));
        tri.push_back(Triangle<1>(3, {0, 4, 5}));
        tri.push_back(Triangle<1>(4, {0, 5, 6}));
        tri.push_back(Triangle<1>(5, {0, 6, 1}));
JasonPries's avatar
JasonPries committed
473

JasonPries's avatar
JasonPries committed
474 475
        inner_region = std::make_shared<Region<2>>(std::vector<size_t>{0, 1, 2, 3, 4, 5});
        inner_boundary = std::make_shared<Boundary<2>>(std::vector<size_t>{1, 2, 3, 4, 5, 6});
JasonPries's avatar
JasonPries committed
476

JasonPries's avatar
JasonPries committed
477 478
        reg.push_back(inner_region);
        bnd.push_back(inner_boundary);
JasonPries's avatar
JasonPries committed
479

480 481
        tri.push_back(Triangle<1>(6, {1, 7, 2}));
        tri.push_back(Triangle<1>(7, {2, 7, 8}));
JasonPries's avatar
JasonPries committed
482

483 484
        tri.push_back(Triangle<1>(8, {2, 8, 3}));
        tri.push_back(Triangle<1>(9, {3, 8, 9}));
JasonPries's avatar
JasonPries committed
485

486 487
        tri.push_back(Triangle<1>(10, {3, 9, 4}));
        tri.push_back(Triangle<1>(11, {4, 9, 10}));
JasonPries's avatar
JasonPries committed
488

489 490
        tri.push_back(Triangle<1>(12, {4, 10, 5}));
        tri.push_back(Triangle<1>(13, {5, 10, 11}));
JasonPries's avatar
JasonPries committed
491

492 493
        tri.push_back(Triangle<1>(14, {5, 11, 6}));
        tri.push_back(Triangle<1>(15, {6, 11, 12}));
JasonPries's avatar
JasonPries committed
494

495 496
        tri.push_back(Triangle<1>(16, {6, 12, 1}));
        tri.push_back(Triangle<1>(17, {1, 12, 7}));
JasonPries's avatar
JasonPries committed
497

JasonPries's avatar
JasonPries committed
498 499
        outer_region = std::make_shared<Region<2>>(std::vector<size_t>{6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17});
        outer_boundary = std::make_shared<Boundary<2>>(std::vector<size_t>{7, 8, 9, 10, 11, 12});
JasonPries's avatar
JasonPries committed
500

JasonPries's avatar
JasonPries committed
501 502
        reg.push_back(outer_region);
        bnd.push_back(outer_boundary);
JasonPries's avatar
JasonPries committed
503 504 505 506 507 508

        femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
    }

    std::vector<XY> node;
    std::vector<Triangle<1>> tri;
JasonPries's avatar
JasonPries committed
509 510
    std::vector<std::shared_ptr<Region<2>>> reg;
    std::vector<std::shared_ptr<Boundary<2>>> bnd;
JasonPries's avatar
JasonPries committed
511

JasonPries's avatar
JasonPries committed
512 513
    std::shared_ptr<Region<2>> inner_region;
    std::shared_ptr<Region<2>> outer_region;
JasonPries's avatar
JasonPries committed
514

JasonPries's avatar
JasonPries committed
515 516
    std::shared_ptr<Boundary<2>> inner_boundary;
    std::shared_ptr<Boundary<2>> outer_boundary;
JasonPries's avatar
JasonPries committed
517

JasonPries's avatar
JasonPries committed
518
    FiniteElementMesh<2, 1> femesh;
JasonPries's avatar
JasonPries committed
519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564

    MaterialProperties Linear1000 = MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1000.0));

    void test_flux_density_element_order_1(Eigen::ArrayXd v, Eigen::ArrayXd Bx, Eigen::ArrayXd By, Eigen::ArrayXd Bmag, Eigen::ArrayXd Bang) {
        double x_q = (1.0 + cos(M_PI / 3.0)) / 3.0;
        double y_q = sin(M_PI / 3.0) / 3.0;
        double r_q = hypot(x_q, y_q);
        double B_q = mu_0 * r_q / 2.0;

        for (size_t i = 0; i != 6; ++i) {
            EXPECT_NEAR(Bmag(i), B_q, FLT_EPSILON);

            double ang = 120.0 + 60.0 * i;
            ang -= 360.0 * (ang > 180.0);

            EXPECT_NEAR(Bang(i), ang, FLT_EPSILON);
        }

        x_q = (1.0 + cos(M_PI / 3.0)) / 2.0;
        y_q = sin(M_PI / 3.0) / 2.0;
        r_q = 2.0 - hypot(x_q, y_q);
        B_q = v(1) / r_q;

        for (size_t i = 6; i < 18; i += 2) {
            EXPECT_NEAR(Bmag(i), B_q, FLT_EPSILON);

            double ang = 120.0 + 60.0 * (i - 6) / 2.0;
            ang -= 360.0 * (ang > 180.0);

            EXPECT_NEAR(Bang(i), ang, FLT_EPSILON);
        }

        x_q = (2.0 * cos(M_PI / 6.0) + 2.0 * cos(M_PI / 6.0 * 3.0)) / 2.0;
        y_q = (2.0 * sin(M_PI / 6.0) + 2.0 * sin(M_PI / 6.0 * 3.0)) / 2.0;
        r_q = hypot(x_q, y_q) - 1.0;
        B_q = v(1) / r_q;

        for (size_t i = 7; i < 18; i += 2) {
            EXPECT_NEAR(Bmag(i), B_q, FLT_EPSILON);

            double ang = 120.0 + 60.0 * (i - 6) / 2.0;
            ang -= 360.0 * (ang > 180.0);

            EXPECT_NEAR(Bang(i), ang, FLT_EPSILON);
        }
    }
JasonPries's avatar
JasonPries committed
565 566
};

JasonPries's avatar
JasonPries committed
567
TEST_F(TwoRegionHexagon, magnetostatic_forcing_air) {
568
    Magnetostatic<2, 1, 1, FieldVariable::A> ms{femesh};
JasonPries's avatar
JasonPries committed
569
    ms.time(0.0);
JasonPries's avatar
JasonPries committed
570 571

    // Test without boundary conditions
JasonPries's avatar
JasonPries committed
572
    ms.build_matrices();
JasonPries's avatar
JasonPries committed
573

JasonPries's avatar
JasonPries committed
574
    auto ff = [](double t) { return 1.0; };
JasonPries's avatar
JasonPries committed
575
    auto f = ms.init_unknown_vector();
JasonPries's avatar
JasonPries committed
576

JasonPries's avatar
JasonPries committed
577 578
    ms.add_current_density(ff, {inner_region});
    ms.add_magnetic_insulation({outer_boundary});
JasonPries's avatar
JasonPries committed
579 580 581 582
    ms.calculate_forcing(f);

    double nodal_current{sqrt(3.0) / 4.0 / 3.0};
    EXPECT_NEAR(f(0), nodal_current * 6.0, FLT_EPSILON);
JasonPries's avatar
JasonPries committed
583 584 585
    for (size_t i = 1; i != 7; ++i) {
        EXPECT_NEAR(f(i), nodal_current * 2.0, FLT_EPSILON);
    }
JasonPries's avatar
JasonPries committed
586

JasonPries's avatar
JasonPries committed
587 588 589
    auto J = ms.init_unknown_matrix();
    auto v = ms.init_unknown_vector();
    auto r = ms.init_unknown_vector();
JasonPries's avatar
JasonPries committed
590 591 592 593 594
    auto Fx = ms.init_element_array();
    auto Fy = ms.init_element_array();
    auto dFxdx = ms.init_element_array();
    auto dFydy = ms.init_element_array();
    auto dFxdy = ms.init_element_array();
JasonPries's avatar
JasonPries committed
595

JasonPries's avatar
JasonPries committed
596 597
    v.setZero();

JasonPries's avatar
JasonPries committed
598
    ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
JasonPries's avatar
JasonPries committed
599

JasonPries's avatar
JasonPries committed
600
    // Test with boundary conditions
JasonPries's avatar
JasonPries committed
601 602
    ms.apply_conditions();

JasonPries's avatar
JasonPries committed
603 604 605 606 607 608
    J = ms.init_unknown_matrix();
    v = ms.init_unknown_vector();
    r = ms.init_unknown_vector();
    f = ms.init_unknown_vector();
    Fx = ms.init_element_vector();
    Fy = ms.init_element_vector();
JasonPries's avatar
JasonPries committed
609 610 611
    dFxdx = ms.init_element_array();
    dFydy = ms.init_element_array();
    dFxdy = ms.init_element_array();
JasonPries's avatar
JasonPries committed
612

JasonPries's avatar
JasonPries committed
613 614 615 616
    v.setZero();

    ms.calculate_forcing(f);
    EXPECT_NEAR(f(0), nodal_current * 6.0, FLT_EPSILON);
JasonPries's avatar
JasonPries committed
617 618 619
    for (size_t i = 1; i != 7; ++i) {
        EXPECT_NEAR(f(i), nodal_current * 2.0, FLT_EPSILON);
    }
JasonPries's avatar
JasonPries committed
620

JasonPries's avatar
JasonPries committed
621
    ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
JasonPries's avatar
JasonPries committed
622

JasonPries's avatar
JasonPries committed
623 624 625 626 627
    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
    ldlt.compute(J);
    ASSERT_EQ(ldlt.info(), Eigen::Success);
    v -= ldlt.solve(r);

JasonPries's avatar
JasonPries committed
628
    // Verify equation is solved
JasonPries's avatar
JasonPries committed
629
    ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
JasonPries's avatar
JasonPries committed
630 631 632 633 634
    for (size_t i = 0; i != r.size(); ++i) {
        EXPECT_NEAR(r(i), 0.0, FLT_EPSILON);
    }

    // Fundamental solution is 0.25*r^2 + C in forced region (0)
JasonPries's avatar
JasonPries committed
635
    EXPECT_NEAR(v(0), v(1) + 0.25 * mu_0, mu_0 * FLT_EPSILON);
JasonPries's avatar
JasonPries committed
636 637 638
    for (size_t i = 1; i != 7; ++i) {
        EXPECT_NEAR(v(1), v(i), FLT_EPSILON);
    }
JasonPries's avatar
JasonPries committed
639

JasonPries's avatar
JasonPries committed
640
    // Test flux-density values
JasonPries's avatar
JasonPries committed
641 642
    Eigen::ArrayXd Bx = ms.derivatives().dy(0).transpose() * v;
    Eigen::ArrayXd By = -ms.derivatives().dx(0).transpose() * v;
JasonPries's avatar
JasonPries committed
643

JasonPries's avatar
JasonPries committed
644 645
    Eigen::ArrayXd Bmag(By.size());
    Eigen::ArrayXd Bang(By.size());
JasonPries's avatar
JasonPries committed
646

JasonPries's avatar
JasonPries committed
647 648 649 650 651
    for (size_t i = 0; i != By.size(); ++i) {
        Bmag(i) = hypot(Bx(i), By(i));
        Bang(i) = atan2(By(i), Bx(i)) * 180.0 / M_PI;
    }

JasonPries's avatar
JasonPries committed
652 653
    test_flux_density_element_order_1(v, Bx, By, Bmag, Bang);
}
JasonPries's avatar
JasonPries committed
654

JasonPries's avatar
JasonPries committed
655
TEST_F(TwoRegionHexagon, magnetostatic_forcing_1000) {
JasonPries's avatar
JasonPries committed
656
    femesh.region(1)->material(Linear1000);
JasonPries's avatar
JasonPries committed
657

658
    Magnetostatic<2, 1, 1, FieldVariable::A> ms{femesh};
JasonPries's avatar
JasonPries committed
659 660
    ms.time(0.0);
    ms.build_matrices();
JasonPries's avatar
JasonPries committed
661

JasonPries's avatar
JasonPries committed
662 663
    auto ff = [](double t) { return 1.0; };

JasonPries's avatar
JasonPries committed
664 665
    ms.add_current_density(ff, {inner_region});
    ms.add_magnetic_insulation({outer_boundary});
JasonPries's avatar
JasonPries committed
666 667 668 669 670 671 672 673 674 675 676 677 678 679

    ms.apply_conditions();

    auto J = ms.init_unknown_matrix();
    auto v = ms.init_unknown_vector();
    auto r = ms.init_unknown_vector();
    auto f = ms.init_unknown_vector();
    auto Fx = ms.init_element_array();
    auto Fy = ms.init_element_array();
    auto dFxdx = ms.init_element_array();
    auto dFydy = ms.init_element_array();
    auto dFxdy = ms.init_element_array();

    v.setZero();
JasonPries's avatar
JasonPries committed
680

JasonPries's avatar
JasonPries committed
681
    ms.calculate_forcing(f);
JasonPries's avatar
JasonPries committed
682

JasonPries's avatar
JasonPries committed
683
    ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
JasonPries's avatar
JasonPries committed
684

JasonPries's avatar
JasonPries committed
685 686 687 688
    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
    ldlt.compute(J);
    ASSERT_EQ(ldlt.info(), Eigen::Success);
    v -= ldlt.solve(r);
JasonPries's avatar
JasonPries committed
689

JasonPries's avatar
JasonPries committed
690 691 692 693
    // Verify equation is solved
    ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
    for (size_t i = 0; i != r.size(); ++i) {
        EXPECT_NEAR(r(i), 0.0, FLT_EPSILON);
JasonPries's avatar
JasonPries committed
694 695
    }

JasonPries's avatar
JasonPries committed
696
    // Fundamental solution is 0.25*r^2 + C in forced region (0)
JasonPries's avatar
JasonPries committed
697
    EXPECT_NEAR(v(0), v(1) + 0.25 * mu_0, mu_0 * FLT_EPSILON);
JasonPries's avatar
JasonPries committed
698 699 700
    for (size_t i = 1; i != 7; ++i) {
        EXPECT_NEAR(v(1), v(i), FLT_EPSILON);
    }
JasonPries's avatar
JasonPries committed
701

JasonPries's avatar
JasonPries committed
702 703 704
    // Test flux-density values
    Eigen::ArrayXd Bx = ms.derivatives().dy(0).transpose() * v;
    Eigen::ArrayXd By = -ms.derivatives().dx(0).transpose() * v;
JasonPries's avatar
JasonPries committed
705

JasonPries's avatar
JasonPries committed
706 707
    Eigen::ArrayXd Bmag(By.size());
    Eigen::ArrayXd Bang(By.size());
JasonPries's avatar
JasonPries committed
708

JasonPries's avatar
JasonPries committed
709 710 711
    for (size_t i = 0; i != By.size(); ++i) {
        Bmag(i) = hypot(Bx(i), By(i));
        Bang(i) = atan2(By(i), Bx(i)) * 180.0 / M_PI;
JasonPries's avatar
JasonPries committed
712
    }
JasonPries's avatar
JasonPries committed
713 714 715 716 717 718 719 720 721 722

    // Flux-density in ur=1000 region should be approximately 1000 times greater
    for (size_t i = 0; i != 6; ++i) {
        for (size_t j = 6; j != 18; ++j) {
            EXPECT_GE(Bmag(j), Bmag(i) * 500);
            EXPECT_LE(Bmag(j), Bmag(i) * 2000);
        }
    }

    test_flux_density_element_order_1(v, Bx, By, Bmag, Bang);
JasonPries's avatar
JasonPries committed
723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739
}

class TwoRegionQuarterOctagon : public ::testing::Test {
public:
    virtual void SetUp() {
        // Inner Region
        node.push_back(XY{0.0, 0.0});

        double r{1.0};
        for (size_t i = 0; i != 3; ++i) {
            double a = (i * M_PI) / 4.0;
            node.push_back(XY{r * std::cos(a), r * std::sin(a)});
        }

        inner_boundary = std::make_shared<Boundary<2>>(std::vector<size_t>{1, 2, 3});
        bnd.push_back(inner_boundary);

740 741
        tri.push_back(Triangle<1>(0, {0, 1, 2}));
        tri.push_back(Triangle<1>(1, {0, 2, 3}));
JasonPries's avatar
JasonPries committed
742 743 744 745 746 747 748 749 750 751 752 753 754 755

        inner_region = std::make_shared<Region<2>>(std::vector<size_t>{0, 1});
        reg.push_back(inner_region);

        // Outer Region
        r = 2.0;
        for (size_t i = 0; i != 3; ++i) {
            double a = (i * M_PI) / 4.0;
            node.push_back(XY{r * std::cos(a), r * std::sin(a)});
        }

        outer_boundary = std::make_shared<Boundary<2>>(std::vector<size_t>{4, 5, 6});
        bnd.push_back(outer_boundary);

756 757 758 759
        tri.push_back(Triangle<1>(2, {1, 4, 2}));
        tri.push_back(Triangle<1>(3, {4, 5, 2}));
        tri.push_back(Triangle<1>(4, {2, 5, 6}));
        tri.push_back(Triangle<1>(5, {2, 6, 3}));
JasonPries's avatar
JasonPries committed
760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788

        outer_region = std::make_shared<Region<2>>(std::vector<size_t>{2, 3, 4, 5});
        reg.push_back(outer_region);

        // Periodic Boundary
        variable_map.push_back(VariableMap{1, 3, 1.0});
        variable_map.push_back(VariableMap{4, 6, 1.0});

        femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
    }

    std::vector<XY> node;
    std::vector<Triangle<1>> tri;
    std::vector<std::shared_ptr<Region<2>>> reg;
    std::vector<std::shared_ptr<Boundary<2>>> bnd;

    std::shared_ptr<Region<2>> inner_region;
    std::shared_ptr<Region<2>> outer_region;

    std::shared_ptr<Boundary<2>> inner_boundary;
    std::shared_ptr<Boundary<2>> outer_boundary;

    std::vector<VariableMap> variable_map;

    FiniteElementMesh<2, 1> femesh;
};

TEST_F(TwoRegionQuarterOctagon, magnetostatic_periodic) {
    // Initialize physics
789
    Magnetostatic<2, 1, 1, FieldVariable::A> ms{femesh};
JasonPries's avatar
JasonPries committed
790 791 792 793 794 795 796 797
    ms.time(0.0);

    // Add conditions
    auto ff = [](double t) { return 1.0; };
    ms.add_current_density(ff, {inner_region});

    ms.add_magnetic_insulation({outer_boundary});

798
    ms.add_periodic_boundary(variable_map, false);
JasonPries's avatar
JasonPries committed
799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857

    ms.build_matrices();

    ms.apply_conditions();

    // Initialize matrices
    auto J = ms.init_unknown_matrix();
    auto v = ms.init_unknown_vector();
    auto r = ms.init_unknown_vector();
    auto f = ms.init_unknown_vector();
    auto Fx = ms.init_element_array();
    auto Fy = ms.init_element_array();
    auto dFxdx = ms.init_element_array();
    auto dFydy = ms.init_element_array();
    auto dFxdy = ms.init_element_array();

    v.setZero();

    ms.calculate_forcing(f);

    double nodal_current{M_SQRT2 / 6.0};
    for (size_t i = 0; i != 3; ++i) {
        EXPECT_NEAR(f(i), nodal_current, FLT_EPSILON);
    }

    ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);

    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
    ldlt.compute(J);
    ASSERT_EQ(ldlt.info(), Eigen::Success);
    v -= ldlt.solve(r);

    // Verify equation is solved
    ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
    for (size_t i = 0; i != r.size(); ++i) {
        EXPECT_NEAR(r(i), 0.0, FLT_EPSILON);
    }

    // Fundamental solution is 0.25*r^2 + C in forced region (0)
    EXPECT_NEAR(v(0), v(1) + 0.25 * mu_0, mu_0 * 0.1);
    for (size_t i = 1; i != 3; ++i) {
        EXPECT_NEAR(v(1), v(i), FLT_EPSILON);
    }

    // Test flux-density values
    Eigen::ArrayXd Bx = ms.derivatives().dy(0).transpose() * v;
    Eigen::ArrayXd By = -ms.derivatives().dx(0).transpose() * v;

    Eigen::ArrayXd Bmag(By.size());
    Eigen::ArrayXd Bang(By.size());

    for (size_t i = 0; i != By.size(); ++i) {
        Bmag(i) = hypot(Bx(i), By(i));
        Bang(i) = atan2(By(i), Bx(i)) * 180.0 / M_PI;
    }

    std::vector<double> b_mag_expected(6, 3.86994e-7);
    std::vector<double> b_ang_expected{112.5, 157.5, 112.5, 112.5, 157.5, 157.5};

JasonPries's avatar
JasonPries committed
858
    for (size_t i = 0; i != 6; ++i) {
JasonPries's avatar
JasonPries committed
859 860 861 862 863 864 865
        EXPECT_NEAR(Bmag[i], b_mag_expected[i], 1e-12);
        EXPECT_NEAR(Bang[i], b_ang_expected[i], 1e-1);
    }
}

TEST_F(TwoRegionQuarterOctagon, magnetostatic_antiperiodic) {
    // Initialize physics
866
    Magnetostatic<2, 1, 1, FieldVariable::A> ms{femesh};
JasonPries's avatar
JasonPries committed
867 868 869 870 871 872 873 874 875
    ms.time(0.0);

    // Add conditions
    auto ff = [](double t) { return 1.0; };
    ms.add_current_density(ff, {inner_region});

    ms.add_magnetic_insulation({outer_boundary});

    variable_map.push_back(VariableMap(0, 0, -1.0));
876
    ms.add_periodic_boundary(variable_map, true);
JasonPries's avatar
JasonPries committed
877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912

    ms.build_matrices();

    ms.apply_conditions();

    // Initialize matrices
    auto J = ms.init_unknown_matrix();
    auto v = ms.init_unknown_vector();
    auto r = ms.init_unknown_vector();
    auto f = ms.init_unknown_vector();
    auto Fx = ms.init_element_array();
    auto Fy = ms.init_element_array();
    auto dFxdx = ms.init_element_array();
    auto dFydy = ms.init_element_array();
    auto dFxdy = ms.init_element_array();

    v.setZero();

    ms.calculate_forcing(f);

    EXPECT_NEAR(f(0), 0.0, FLT_EPSILON);
    EXPECT_NEAR(f(1), M_SQRT2 / 6.0, FLT_EPSILON);

    ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);

    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
    ldlt.compute(J);
    ASSERT_EQ(ldlt.info(), Eigen::Success);
    v -= ldlt.solve(r);

    // Verify equation is solved
    ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
    for (size_t i = 0; i != r.size(); ++i) {
        EXPECT_NEAR(r(i), 0.0, FLT_EPSILON);
    }

JasonPries's avatar
JasonPries committed
913
    std::vector<double> v_expected{0.0, 6.60365e-8};
JasonPries's avatar
JasonPries committed
914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933

    for (size_t i = 0; i != 2; ++i) {
        EXPECT_NEAR(v[i], v_expected[i], 1e-13);
    }

    // Test flux-density values
    Eigen::ArrayXd Bx = ms.derivatives().dy(0).transpose() * v;
    Eigen::ArrayXd By = -ms.derivatives().dx(0).transpose() * v;

    Eigen::ArrayXd Bmag(By.size());
    Eigen::ArrayXd Bang(By.size());

    for (size_t i = 0; i != By.size(); ++i) {
        Bmag(i) = hypot(Bx(i), By(i));
        Bang(i) = atan2(By(i), Bx(i)) * 180.0 / M_PI;
    }

    std::vector<double> b_mag_expected{9.33897e-8, 9.33897e-8, 9.33897e-8, 7.14774e-8, 7.14774e-8, 9.33897e-8};
    std::vector<double> b_ang_expected{0.0, -90, 0.0, 112.5, 157.5, -90.0};

JasonPries's avatar
JasonPries committed
934
    for (size_t i = 0; i != 6; ++i) {
JasonPries's avatar
JasonPries committed
935 936 937
        EXPECT_NEAR(Bmag[i], b_mag_expected[i], 1e-13);
        EXPECT_NEAR(Bang[i], b_ang_expected[i], FLT_EPSILON);
    }
JasonPries's avatar
JasonPries committed
938 939 940 941 942 943 944 945 946 947 948 949 950 951
}

class TwoRegionHexagonRotation : public ::testing::Test {
public:
    virtual void SetUp() {
        // Inner Region
        node.push_back(XY{0.0, 0.0});

        double r{1.0};
        for (size_t i = 0; i != 6; ++i) {
            double a = (i * M_PI) / 3.0;
            node.push_back(XY{r * std::cos(a), r * std::sin(a)});
        }

952 953 954 955 956 957
        tri.push_back(Triangle<1>(0, {0, 1, 2}));
        tri.push_back(Triangle<1>(1, {0, 2, 3}));
        tri.push_back(Triangle<1>(2, {0, 3, 4}));
        tri.push_back(Triangle<1>(3, {0, 4, 5}));
        tri.push_back(Triangle<1>(4, {0, 5, 6}));
        tri.push_back(Triangle<1>(5, {0, 6, 1}));
JasonPries's avatar
JasonPries committed
958

JasonPries's avatar
JasonPries committed
959 960
        region_ip = std::make_shared<Region<2>>(std::vector<size_t>{0, 1, 2});
        region_in = std::make_shared<Region<2>>(std::vector<size_t>{3, 4, 5});
JasonPries's avatar
JasonPries committed
961 962
        boundary_io = std::make_shared<Boundary<2>>(std::vector<size_t>{1, 2, 3, 4, 5, 6});

JasonPries's avatar
JasonPries committed
963 964
        reg.push_back(region_ip);
        reg.push_back(region_in);
JasonPries's avatar
JasonPries committed
965 966 967 968 969 970 971 972 973 974 975 976 977 978 979
        bnd.push_back(boundary_io);

        // Outer Region
        r = 1.0;
        for (size_t i = 0; i != 6; ++i) {
            double a = (i * M_PI) / 3.0;
            node.push_back(XY{r * std::cos(a), r * std::sin(a)});
        }

        r = 2.0;
        for (size_t i = 0; i != 6; ++i) {
            double a = (i * M_PI) / 3.0 + M_PI / 6.0;
            node.push_back(XY{r * std::cos(a), r * std::sin(a)});
        }

980 981
        tri.push_back(Triangle<1>(6, {7, 13, 8}));
        tri.push_back(Triangle<1>(7, {8, 13, 14}));
JasonPries's avatar
JasonPries committed
982

983 984
        tri.push_back(Triangle<1>(8, {8, 14, 9}));
        tri.push_back(Triangle<1>(9, {9, 14, 15}));
JasonPries's avatar
JasonPries committed
985

986 987
        tri.push_back(Triangle<1>(10, {9, 15, 10}));
        tri.push_back(Triangle<1>(11, {10, 15, 16}));
JasonPries's avatar
JasonPries committed
988

989 990
        tri.push_back(Triangle<1>(12, {10, 16, 11}));
        tri.push_back(Triangle<1>(13, {11, 16, 17}));
JasonPries's avatar
JasonPries committed
991

992 993
        tri.push_back(Triangle<1>(14, {11, 17, 12}));
        tri.push_back(Triangle<1>(15, {12, 17, 18}));
JasonPries's avatar
JasonPries committed
994

995 996
        tri.push_back(Triangle<1>(16, {12, 18, 7}));
        tri.push_back(Triangle<1>(17, {7, 18, 13}));
JasonPries's avatar
JasonPries committed
997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013

        region_o = std::make_shared<Region<2>>(std::vector<size_t>{6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17});
        boundary_oi = std::make_shared<Boundary<2>>(std::vector<size_t>{7, 8, 9, 10, 11, 12});
        boundary_oo = std::make_shared<Boundary<2>>(std::vector<size_t>{13, 14, 15, 16, 17, 18});

        reg.push_back(region_o);
        bnd.push_back(boundary_oi);
        bnd.push_back(boundary_oo);

        femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
    }

    std::vector<XY> node;
    std::vector<Triangle<1>> tri;
    std::vector<std::shared_ptr<Region<2>>> reg;
    std::vector<std::shared_ptr<Boundary<2>>> bnd;

JasonPries's avatar
JasonPries committed
1014 1015
    std::shared_ptr<Region<2>> region_ip;
    std::shared_ptr<Region<2>> region_in;
JasonPries's avatar
JasonPries committed
1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083
    std::shared_ptr<Region<2>> region_o;

    std::shared_ptr<Boundary<2>> boundary_io;
    std::shared_ptr<Boundary<2>> boundary_oi;
    std::shared_ptr<Boundary<2>> boundary_oo;

    FiniteElementMesh<2, 1> femesh;

    MaterialProperties Linear1000 = MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1000.0));

    void test_flux_density_element_order_1(Eigen::ArrayXd v, Eigen::ArrayXd Bx, Eigen::ArrayXd By, Eigen::ArrayXd Bmag, Eigen::ArrayXd Bang) {
        double x_q = (1.0 + cos(M_PI / 3.0)) / 3.0;
        double y_q = sin(M_PI / 3.0) / 3.0;
        double r_q = hypot(x_q, y_q);
        double B_q = mu_0 * r_q / 2.0;

        for (size_t i = 0; i != 6; ++i) {
            EXPECT_NEAR(Bmag(i), B_q, FLT_EPSILON);

            double ang = 120.0 + 60.0 * i;
            ang -= 360.0 * (ang > 180.0);

            if (std::abs(180.0 + Bang(i)) < 0.1) {
                Bang(i) = -Bang(i);
            }

            EXPECT_NEAR(Bang(i), ang, FLT_EPSILON);
        }

        x_q = (1.0 + cos(M_PI / 3.0)) / 2.0;
        y_q = sin(M_PI / 3.0) / 2.0;
        r_q = 2.0 - hypot(x_q, y_q);
        B_q = v(1) / r_q;

        for (size_t i = 6; i < 18; i += 2) {
            EXPECT_NEAR(Bmag(i), B_q, FLT_EPSILON);

            double ang = 120.0 + 60.0 * (i - 6) / 2.0;
            ang -= 360.0 * (ang > 180.0);

            if (std::abs(180.0 + Bang(i)) < 0.1) {
                Bang(i) = -Bang(i);
            }

            EXPECT_NEAR(Bang(i), ang, FLT_EPSILON);
        }

        x_q = (2.0 * cos(M_PI / 6.0) + 2.0 * cos(M_PI / 6.0 * 3.0)) / 2.0;
        y_q = (2.0 * sin(M_PI / 6.0) + 2.0 * sin(M_PI / 6.0 * 3.0)) / 2.0;
        r_q = hypot(x_q, y_q) - 1.0;
        B_q = v(1) / r_q;

        for (size_t i = 7; i < 18; i += 2) {
            EXPECT_NEAR(Bmag(i), B_q, FLT_EPSILON);

            double ang = 120.0 + 60.0 * (i - 6) / 2.0;
            ang -= 360.0 * (ang > 180.0);

            if (std::abs(180.0 + Bang(i)) < 0.1) {
                Bang(i) = -Bang(i);
            }

            EXPECT_NEAR(Bang(i), ang, FLT_EPSILON);
        }
    }
};

TEST_F(TwoRegionHexagonRotation, magnetostatic_forcing) {
JasonPries's avatar
JasonPries committed
1084
    // Tests sliding interface condition for implementing virtual rotational motion
1085
    Magnetostatic<2, 1, 1, FieldVariable::A> ms{femesh};
JasonPries's avatar
JasonPries committed
1086